arXiv:physics/0008016  vl  7  Aug  2000 


Heart  Rate  Variability:  Measures  and  Models 


Malvin  C.  Tcich,  Steven  B.  Lowen,  Bradley  M.  Jost,  and  Karin  Vibe-Rhcymer 
Department  of  Electrical  and  Computer  Engineering,  Boston  University 
8  Saint  Mary’s  Street,  Boston,  MA  02215 
and 

Conor  Heneghan 

Department  of  Electronic  and  Electrical  Engineering, 

University  College  Dublin,  Belfield,  Dublin  4,  Ireland 


1 


Report  Documentation  Page 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 

VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 

1.  REPORT  DATE 

07  JAN  2005 

2.  REPORT  TYPE 

N/A 

3.  DATES  COVERED 

4.  TITLE  AND  SUBTITLE 

5a.  CONTRACT  NUMBER 

Heart  Rate  Variability:  Measures  and  Models 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Department  of  Electrical  and  Computer  Engineering,  Boston  University  8 
Saint  Mary’s  Street,  Boston,  MA  02215;  Department  of  Electronic  and 
Electrical  Engineering,  University  College  Dublin,  Bel  eld,  Dublin  4, 
Ireland 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release,  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

See  also  ADM001750,  Wavelets  and  Multifractal  Analysis  (WAMA)  Workshop  held  on  19-31  July  2004., 

The  original  document  contains  color  images. 

14.  ABSTRACT 

15.  SUBJECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 
ABSTRACT 

uu 

18.  NUMBER 
OF  PAGES 

84 

19a.  NAME  OF 
RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


Contents 


1  Introduction 


2  Methods  and  Measures 


2.1  The  Heartbeat  Sequence  as  a  Point  Process 


2.1.1  Conventional  Point  Processes 


2.1.2  Fractal  and  Fractal-Rate  Point  Processes 


2.2  Standard  Frequency-Domain  Measures 


2.3  Standard  Time-Domain  Measures 


IA  Other  Standard  Measures 


2.5  Novel  Scalc-Dcpcndcnt  Measures! 


2.5.1  Allan  Factor  \A[T] 


2.5.2  Wavelet-Transform  Standard  Deviation  Icr 


[m 


2.5.3  Relationship  of  Wavelet 

O-wav(m) 

and  Spectral  Measures 

MR. 

2.5.4  Detrended  Fluctuation  Analysis  |DFA(m) 


2.6  Scale-Independent  Measures 


2.6.1  Detrended-Fluctuation-Analysis  Power-Law  Exponent  (ap) 


2.6.2  Wavelet-Transform  Power-Law  Exponent  [aw) 

2.6.3  Periodogram  Power-Law  Exponent  (og)  .  .  . 


2.6.4  Allan-Factor  Power-Law  Exponent,  (oa) 


2.6.5  Rescaled-Range-Analysis  Power-Law  Exponent  [an) 


2.7  Estimating  the  Performance  of  a  Measure  .  . 
27T.1  Statistical  Significance:  p,  d' ,  h,  and  a 


5 

6 

6 

7 

9 

10 

12 

12 

13 

13 

14 
16 

17 

18 
18 
18 
18 
19 

19 

20 
20 


2 


21 


2.7.2  Positive  and  Negative  Predictive  Values 


2.7.3  Receiver-Operating-Characteristic  (ROC)  Analysis 


3  Discriminating  Heart-Failure  Patients  from  Normal  Subjects 


3.5.2  ROC-Area  Curves 


3.6  Comparing  the  Measures  for  Various  Data  Lengths  .  . 
3U.1  Scale-independent  vs  Scale-Dependent  Measures 


3.6.2  Computation  Times  of  the  Various  Measures 

3.6.3  Comparing  the  Most  Effective  Measures  .  . 


3.6.4  Selecting  the  Best  Measure^ 


4  Markers  for  Other  Cardiac  Pathologies 


5  Does  Deterministic  Chaos  Play  a  Role  in  Heart  Rate  Variability? 

5.1  Methods . 


5.1.1  Phase-Space  Reconstruction 

5.1.2  Removing  Correlations  in  the  Data 

5.1.3  Surrogate  Data  Analysis 

5.2  Absence  of  Chaos 


21 

22 

22 

23 

24 

24 

26 

26 

27 

28 

29 

30 

30 

32 

33 

34 

35 

35 

36 

36 

37 


3 


6  Mathematical  Models  for  Heart  Rate  Variability 


|6.1  Integrate-and-Fire  Model . 

Kernel  ol  the  intcgratc-and-Pirc  Model 


6.3.1  Simulating  the  Jittered  Integrate-and-Fire  Point.  Process 


6.3.2  Statistics  ol  the  Simulated  Point  Process  tor  Normal  Subjects  and  (JHF 

Patients 

6.3.3  Simulated  individual  Value  Plots  and  ROC- Area  Curves 


6.3.4  Limitations  of  the  Jittered  Integrate-and-Fire  Model 
3.4  Toward  an  Improved  Model  ol  Heart  Kate  Variability  .  .  . 


38 

39 

39 

39 

40 

41 

41 

42 

42 

43 

44 

45 

46 

48 

49 

49 

53 

56 


4 


1  Introduction 


The  human  heart  generates  the  quintessential  biological  signal:  the  heartbeat.  A  recording 
of  the  cardiac-induced  skin  potentials  at  the  body’s  surface,  an  electrocardiogram  (ECG), 
reveals  information  about  atrial  and  ventricular  electrical  activity.  Abnormalities  in  the 
temporal  durations  of  the  segments  between  deflections,  or  of  the  intervals  between  waves  in 
the  ECG,  as  well  as  their  relative  heights,  serve  to  expose  and  distinguish  cardiac  dysfunc¬ 
tion.  Because  the  electrical  activity  of  the  human  heart  is  influenced  by  many  physiological 
mechanisms,  electrocardiography  has  become  an  invaluable  tool  for  the  diagnosis  of  a  variety 
of  pathologies  that  affect  the  cardiovascular  system  [jTJ.  Electrocardiologists  have  come  to 
excel  at  visually  interpreting  the  detailed  form  of  the  ECG  wave  pattern  and  have  become 
adept  at  differential  diagnoses. 

Readily  recognizable  features  of  the  ECG  wave  pattern  are  designated  by  the  letters  P- 
QRS-T;  the  wave  itself  is  often  referred  to  as  the  QRS  complex.  Aside  from  the  significance 
of  various  features  of  the  QRS  complex,  the  timing  of  the  sequence  of  QRS  complexes  over 
tens,  hundreds,  and  thousands  of  heartbeats  is  also  significant.  These  inter-complex  times 
are  readily  measured  by  recording  the  occurrences  of  the  peaks  of  the  large  R  waves  which 
are,  perhaps,  the  most  distinctive  feature  of  the  normal  ECG. 

In  this  Chapter  we  focus  on  various  measures  of  the  fluctuations  of  this  sequence  of 
interbeat  intervals  and  how  such  fluctuations  can  be  used  to  assess  the  presence  or  likelihood 
of  cardiovascular  disease  [|j.  This  approach  has  come  to  be  called  heart  rate  variability 
(HRV)  analysis  ||  [|  even  when  it  is  the  time  intervals  whose  fluctuations  are  studied 
(heart  rate  has  units  of  inverse  time  rather  than  time).  HRV  analysis  serves  as  a  marker  for 
cardiovascular  disease  because  cardiac  dysfunction  is  often  manifested  by  systematic  changes 
in  the  variability  of  the  RR-interval  sequence  relative  to  that  of  normal  controls  [0.  SSI-  A 
whole  host  of  HRV  measures,  some  scale-dependent  and  others  scale-independent,  have  been 
developed  and  examined  over  the  years  in  an  effort  to  develop  readily  available,  inexpensive, 
and  noninvasive  measures  of  cardiovascular  function. 

We  examine  sixteen  HRV  measures  and  their  suitability  for  correctly  classifying  ECG 
records  of  various  lengths  as  normal  or  revealing  the  presence  of  cardiac  dysfunction.  Par¬ 
ticular  attention  is  devoted  to  HRV  measures  that  are  useful  for  discriminating  congestive- 
heart-failure  patients  from  normal  subjects.  Using  receiver-operating-characteristic  (ROC) 
analysis  we  demonstrate  that  scale-dependent  HRV  measures  (e.g.,  wavelet  and  spectral 
measures)  are  substantially  superior  to  scale-independent  measures  (such  as  wavelet  and 
spectral  fractal  exponents)  for  discriminating  these  two  classes  of  data  over  a  broad  range 
of  record  lengths.  The  wavelet-transform  standard  deviation  at  a  scale  near  32  heartbeat 
intervals,  and  its  spectral  counterpart  near  1/32  cycles/interval,  turn  out  to  provide  reliable 
results  using  ECG  records  just  minutes  long. 
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A  long-standing  issue  of  importance  in  cardiac  physiology  is  the  determination  of  whether 
the  normal  RR  sequence  arises  from  a  chaotic  attractor  or  has  an  underlying  stochastic  origin 
[6;].  We  present  a  phase-space  analysis  in  which  differences  between  adjacent  RR  intervals 
are  embedded.  This  has  the  salutary  effect  of  removing  most  of  the  correlation  in  the  time 
series,  which  is  well-known  to  be  deleterious  to  the  detection  of  underlying  deterministic 
dynamics.  We  demonstrate  that  RR  sequences,  from  normal  subjects  and  from  patients 
with  cardiac  dysfunction  alike,  have  stochastic  rather  than  deterministic  origins,  in  accord 
with  our  earlier  conclusions  [7:.  . 

Finally  we  develop  a  mathematical  point  process  that  emulates  the  human  heartbeat 
time  series  for  both  normal  subjects  and  heart-failure  patients.  Using  simulations,  we  show 
that  a  jittered  integrate-and-hre  model  built  around  a  fractal-Gaussian-noise  kernel  provides 
a  realistic,  though  not  perfect,  simulation  of  real  heartbeat  sequences.  A  construct  of  this 
kind  may  well  be  useful  in  a  number  of  venues,  including  pacemaker  excitation. 


2  Methods  and  Measures 

2.1  The  Heartbeat  Sequence  as  a  Point  Process 

The  statistical  behavior  of  the  sequence  of  heartbeats  can  be  studied  by  replacing  the  complex 
waveform  of  an  individual  heartbeat  recorded  in  the  ECG  (an  entire  QRS-complex)  with  the 
time  of  occurrence  of  the  contraction  (the  time  of  the  peak  of  the  R  phase),  which  is  a 
single  number  [|],  ||] .  In  mathematical  terms,  the  heartbeat  sequence  is  then  modeled  as  an 
unmarked  point  process.  This  simplification  greatly  reduces  the  computational  complexity 
of  the  problem  and  permits  us  to  use  the  substantial  methodology  that  exists  for  point 
processes  13  0.0. 

The  occurrence  of  a  contraction  at  time  U  is  therefore  simply  represented  by  an  impulse 
S(t  —  tf)  at  that  time,  where  5  is  the  Dirac  delta  function,  so  that  the  sequence  of  heartbeats 
is  represented  by 

Ht)  =  ^2s(t-u).  (i) 

i 

A  realization  of  a  point  process  is  specified  by  the  set  of  occurrence  times  {U}  of  the  events. 
A  single  realization  of  the  data  is  often  all  that  is  available  to  the  observer  so  that  the 
identification  of  the  point  process,  and  the  elucidation  of  the  mechanisms  that  underlie  it, 
must  be  gleaned  from  this  one  realization. 

One  way  in  which  the  information  in  an  experimental  point  process  can  be  made  more 
digestible  is  to  reduce  the  data  into  a  statistic  that  emphasizes  a  particular  aspect  of  the 
data  (at  the  expense  of  other  features).  These  statistics  fall  into  two  broad  classes  which 
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derive  from  the  sequence  of  interevent  intervals  and  the  sequence  of  counts,  as  illustrated  in 
Fig.  10,0. 

Figure  1  illustrates  how  an  electrocardiogram  may  be  analyzed  to  obtain  the  sequence  of 
interbeat  intervals  as  well  as  the  sequence  of  counts.  Fig.  1(a)  illustrates  an  ECG  (sequence 
of  QRS  complexes)  recorded  from  a  patient.  The  R  waves  are  schematically  represented  by 
a  sequence  of  vertical  lines,  as  shown  in  Fig.  1(b).  The  time  between  the  first  two  R  waves 
is  Ti,  the  first  RR  (or  interbeat)  interval,  as  indicated  by  the  horizontal  arrows  in  this  figure. 
The  time  between  the  second  and  third  R  waves  is  72,  and  so  forth.  In  Fig.  1(c),  the  time 
axis  is  divided  into  equally  spaced,  contiguous  time  windows,  each  of  duration  T  seconds, 
and  the  (integer)  number  of  R  waves  that  fall  in  the  ith  window  is  counted  and  denoted 
Ni.  This  sequence  {Ni}  forms  a  discrete-time  random  counting  process  of  nonnegative 
integers.  Varying  the  duration  T  yields  a  family  of  sequences  {lVj}(T).  The  RR  intervals 
{rj}  themselves  also  form  a  sequence  of  positive  real-valued  random  numbers,  which  is  shown 
schematically  in  Fig.  1(d).  Here  the  abscissa  is  the  interval  number,  which  is  not  a  simple 
function  of  time. 

In  this  section  we  examine  several  statistical  measures  (including  some  that  are  novel)  to 
characterize  these  stochastic  processes;  the  development  is  assisted  by  an  understanding  of 
point  processes. 


2.1.1  Conventional  Point  Processes 


The  homogeneous  Poisson  point  process,  perhaps  the  simplest  of  all  stochastic  point  pro¬ 
cesses,  is  described  by  a  single  parameter,  the  rate  A.  This  point  process  is  memoryless:  the 
occurrence  of  an  event  at  any  time  to  is  independent  of  the  presence  (or  absence)  of  events 
at  other  times  t  7^  t0.  Because  of  this  property,  both  the  intervals  {rj}  and  counts  {Ni} 
form  sequences  of  independent,  identically  distributed  random  variables.  The  homogeneous 
Poisson  point  process  is  therefore  completely  characterized  by  the  interevent-interval  distri¬ 
bution  (also  referred  to  as  the  interbeat-interval  histogram),  which  is  exponential,  or  the 
event-number  distribution  (also  referred  to  as  the  counting  distribution),  which  is  Poisson, 
together  with  the  property  of  being  independent.  This  process  serves  as  a  benchmark  against 
which  other  point  processes  are  measured;  it  therefore  plays  the  role  that  the  white  Gaussian 
process  enjoys  in  the  realm  of  continuous-time  stochastic  processes. 


A  related  point  process  is  the  nonparalyzable  fixed-dead-time-modified  Poisson  point 
process,  a  close  cousin  of  the  homogeneous  Poisson  point  process  that  differs  only  by  the 
imposition  of  a  dead-time  (refractory)  interval  after  the  occurrence  of  each  event,  during 
which  other  events  are  prohibited  from  occurring  0.0.  Another  cousin  is  the  gamma-r 
renewal  process  which,  for  integer  r,  is  generated  from  an  homogeneous  Poisson  point  process 
by  permitting  every  rth  event  to  survive  while  deleting  all  intermediate  events  [ITUl  [HI .  Both 
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the  dead-time-modified  Poisson  point  process  and  the  gamma-r  renewal  process  require  two 
parameters  for  their  description. 


Some  point  processes  exhibit  no  dependencies  among  their  interevent  intervals  at  the 
outset,  in  which  case  the  sequence  of  interevent  intervals  forms  a  sequence  of  identically 
distributed  random  variables  and  the  point  process  is  completely  specified  by  its  interevent- 
interval  histogram,  i.e. ,  its  first-order  statistic.  Such  a  process  is  called  a  renewal  process 
T0|,  a  definition  motivated  by  the  replacement  of  failed  parts,  each  replacement  of  which 
forms  a  renewal  of  the  point  process.  Both  examples  of  point  processes  presented  above 
belong  to  the  class  of  renewal  point  processes. 


The  interevent-interval  histogram  is,  perhaps,  the  most  commonly  used  of  all  statistical 
measures  of  point  processes  in  the  life  sciences.  The  interevent-interval  histogram  estimates 
the  interevent-interval  probability  density  function  pT  (r)  by  computing  the  relative  frequency 
of  occurrence  of  interevent  intervals  as  a  function  of  interval  size.  Its  construction  involves  the 
loss  of  interval  ordering,  and  therefore  of  information  about  dependencies  among  intervals; 
a  reordering  of  the  sequence  does  not  alter  the  interevent-interval  histogram  since  the  order 
plays  no  role  in  the  relative  frequency  of  occurrence. 


The  interevent-interval  probability  density  function  for  the  homogeneous  Poisson  point 
process  assumes  the  exponential  form 


pT(r)  =  Aexp(— Ar) 


(2) 


where  A  is  the  mean  number  of  events  per  unit  time.  The  interevent-interval  mean  and 
variance  are  readily  calculated  to  be  E[r]  =  /0°°  rpT(r)dr  =  1/A  and  Var(r)  =  E[r2]  — E2[r]  = 
1/A2,  respectively,  where  E[-]  represents  expectation  over  the  quantity  inside  the  brackets. 
The  interevent-interval  probability  density  function  for  the  dead-time-modified  Poisson  point 
process  exhibits  the  same  exponential  form  as  for  the  homogeneous  Poisson  point  process, 
but  is  truncated  at  short  interevent  intervals  as  a  result  of  the  dead  time 
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Pr(T|“{  Aexp[— A(t  —  Td)] 


T  <Td 
T>Td 


(3) 


Here  rd  is  the  dead  time  and  A  is  the  rate  of  the  process  before  dead  time  is  imposed. 


If  a  process  is  nonrenewal,  so  that  dependencies  exist  among  its  interevent  intervals, 
then  the  interevent-interval  histogram  does  not  completely  characterize  the  process  [O .  In 


this  case,  measures  that  reveal  the  nature  of  the  dependencies  provide  information  that  is 
complementary  to  that  contained  in  the  interevent-interval  histogram.  The  heartbeat  time 
series  is  such  a  nonrenewal  process. 


2.1.2  Fractal  and  Fractal-Rate  Point  Processes 


The  complete  characterization  of  a  stochastic  process  involves  a  description  of  all  possible 
joint  probabilities  of  the  various  events  occurring  in  the  process.  Different  statistics  provide 
complementary  views  of  the  process;  no  single  statistic  can  in  general  describe  a  stochastic 
process  completely.  Fractal  stochastic  processes  exhibit  scaling  in  their  statistics.  Such 
scaling  leads  naturally  to  power-law  behavior,  as  demonstrated  in  the  following.  Consider  a 
statistic  w,  such  as  the  Allan  factor  for  long  counting  times  (see  Sec.  US),  which  depends 
continuously  on  the  scale  x  over  which  measurements  are  taken  [[IT].  |T7|] .  Suppose  changing 
the  scale  by  any  factor  a  effectively  scales  the  statistic  by  some  other  factor  g(a),  related  to 
the  factor  but  independent  of  the  original  scale: 


w(ax)  =  g(a)w(x) 


(4) 


The  only  nontrivial  solution  of  this  scaling  equation,  for  real  functions  and  arguments,  that 
is  independent  of  a  and  x  is 


w(x)  =  bg{x)  with  g{x)  =  xc 


(5) 


for  some  constants  b  and  c  m  B0-  Thus  statistics  with  power-law  forms  are  closely 


related  to  the  concept  of  a  fractal  ||T9|,  [20].  pljl .  The  particular  case  of  fixed  a  admits  a  more 
general  solution  [P2|: 

g(x;  a)  =  xc  cos[27t  ln(x) /  ln(a)] .  (6) 


Consider  once  again,  for  example,  the  interevent-interval  histogram.  This  statistic  high¬ 
lights  the  behavior  of  the  times  between  adjacent  events,  but  reveals  none  of  the  information 
contained  in  the  relationships  among  these  times,  such  as  correlation  between  adjacent  time 
intervals.  If  the  interevent-interval  probability  density  function  follows  the  form  of  Eq.  (p|) 
so  that  p(r)  ~  rc  over  a  certain  range  of  r  where  c  <  —  1,  the  process  is  known  as  a  fractal 
renewal  point  process  0  0.  a  form  of  fractal  stochastic  process. 


A  number  of  statistics  may  be  used  to  describe  a  fractal  stochastic  point  process,  and 
each  statistic  which  scales  will  in  general  have  a  different  scaling  exponent  c.  Each  of  these 
exponents  can  be  simply  related  to  a  more  general  parameter  a,  the  fractal  exponent,  where 
the  exact  relation  between  these  two  exponents  will  depend  upon  the  statistic  in  question. 
For  example,  the  exponent  c  of  the  interevent-interval  probability  density  function  defined 
above  is  related  to  the  fractal  exponent  a  by  c  =  —  (1  +  a).  As  the  fractal  exponent  is  a 
constant  that  describes  the  overall  scaling  behavior  of  a  statistic,  it  does  not  depend  on  the 
particular  scale  and  is  therefore  scale  independent.  Scale-independent  measures  are  discussed 
in  subsections  and  [ht>l[ 


Sample  functions  of  the  fractal  renewal  point  process  are  true  fractals;  the  expected 
value  of  their  generalized  dimensions  assumes  a  nonintegral  value  between  the  topological 
dimension  (zero)  and  the  Euclidean  dimension  (unity)  0- 
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The  sequence  of  unitary  events  observed  in  many  biological  and  physical  systems,  such  as 
the  heartbeat  sequence,  do  not  exhibit  power-law-distributed  interevent-interval  histograms 
but  nevertheless  exhibit  scaling  in  other  statistics.  These  processes  therefore  have  integral 
generalized  dimensions  and  are  consequently  not  true  fractals.  They  may  nevertheless  be 
endowed  with  rate  functions  that  are  either  fractals  or  their  increments:  fractal  Brownian 
motion,  fractal  Gaussian  noise,  or  other  related  processes.  Therefore,  such  point  processes 
are  more  properly  termed  fractal-rate  stochastic  point  processes  [Q-  It  can  be  shown  by 
surrogate  data  methods,  e.g.,  shuffling  the  order  of  the  intervals  (see  Sec.  |5.1.3|),  that  it  is 
the  ordering  and  not  the  relative  interval  sizes  that  distinguish  these  point  processes  [0 . 


2.2  Standard  Frequency-Domain  Measures 

A  number  of  HRV  measures  have  been  used  as  standards  in  cardiology,  both  for  purposes  of 
physiological  interpretation  and  for  clinical  diagnostic  applications  [3;].  We  briefly  describe 
some  of  the  more  commonly  used  measures  that  we  include  in  this  chapter  for  comparison 
with  several  novel  measures  that  have  been  recently  developed. 

Fourier  transform  techniques  provide  a  method  for  quantifying  the  correlation  properties 
of  a  stochastic  process  through  spectral  analysis.  Two  definitions  of  power  spectral  density 
have  been  used  in  the  analysis  of  HRV  p}.  A  rate-based  power  spectral  density  S\(f)  is 
obtained  by  deriving  an  underlying  random  continuous  process  A (t),  the  heart  rate,  based 
on  a  transformation  of  the  observed  RR  interbeat  intervals.  The  power  spectral  density  of 
this  random  process  is  well  defined,  and  standard  techniques  may  be  used  for  estimating  the 
power  spectral  density  from  a  single  observation  of  A (t).  An  advantage  of  this  technique  is 
that  the  power  spectral  density  thus  calculated  has  temporal  frequency  as  the  independent 
variable,  so  that  spectral  components  can  be  interpreted  in  terms  of  underlying  physiological 
processes  with  known  timescales.  The  power  spectral  density  itself  is  usually  expressed  in 
units  of  sec-1.  However,  the  choice  of  how  to  calculate  A (t),  the  underlying  rate  function, 
may  influence  the  calculated  power  spectral  density. 

The  second  spectral  measure  that  is  more  widely  used  is  an  interval- based  power  spec¬ 
tral  density  ST(f )  that  is  directly  calculated  from  measured  RR  interbeat  intervals  without 
transformation  J9|.  In  this  case  the  intervals  are  treated  as  discrete-index  samples  of  an 
underlying  random  process,  and  there  is  no  intermediate  calculation  of  an  underlying  rate 
function.  The  power  spectral  density  in  this  case  has  cycles/interval  as  the  independent 
variable,  and  therefore  has  units  of  sec2/interval. 

The  two  types  of  power  spectral  densities  are  easily  confused  and  care  must  be  taken  in 
their  interpretation.  For  example,  one  could  mistakenly  interpret  the  abscissa  of  an  interval- 
based  power  spectral  density  plot  as  being  equivalent  to  temporal  frequency  (e.g.,  cycles/sec). 
While  this  is  generally  incorrect,  for  point  processes  whose  interevent-interval  coefficient  of 
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variation  is  relatively  small  ||,  the  interval-based  and  rate-based  power  spectral  density 
plots  can  be  made  approximately  equivalent  by  converting  the  interval-based  frequency  /int 
(in  cycles/interval)  to  the  time-based  frequency  ft ime  (in  cycles/sec)  using 

/time  =  /int/E[r].  (7) 

For  typical  interbeat-interval  sequences,  the  coefficient  of  variation  is  indeed  relatively  small 
and  this  conversion  can  be  carried  out  without  the  introduction  of  significant  error  In 
the  remainder  of  this  Chapter  we  work  principally  with  the  interval-based  power-spectral 
density.  We  use  the  notation  /  =  /int  for  the  interval-based  frequency  (cycles/interval)  and 
retain  the  notation  ft ime  for  temporal  frequency  (cycles/sec). 

We  make  use  of  a  non-parametric  technique  for  estimating  the  spectral  density.  A  simple 
reliable  method  for  estimating  the  power  spectral  density  of  a  process  from  a  set  of  discrete 
samples  {rj}  is  to  calculate  the  averaged  periodogram  ^4],  EH-  The  data  is  first  divided 
into  K  non-overlapping  blocks  of  L  samples.  After  the  optional  use  of  a  Hanning  window, 
the  discrete  Fourier  transform  of  each  block  is  calculated  and  squared.  The  results  are  then 
averaged  to  form  the  estimate 


Sr(f)^fzm)  r.  (s) 

k= 1 

Here  fk{f)  is  the  discrete  Fourier  transform  of  the  7th  block  of  data  and  the  hat  explic¬ 
itly  indicates  that  we  are  dealing  with  an  estimate  of  ST(f),  which  is  called  an  averaged 
periodogram. 


The  periodogram  covers  a  broad  range  of  frequencies  which  can  be  divided  into  bands 
that  are  relevant  to  the  presence  of  various  cardiac  pathologies.  The  power  within  a  band 
is  calculated  by  integrating  the  power  spectral  density  over  the  associated  frequency  range. 
Some  commonly  used  measures  in  HRV  are  [|J : 


VLF.  The  power  in  the  very-low-frequency  range:  0.003-0.04  cyclcs/interval.  Physiological 
correlates  of  the  VLF  band  have  not  been  specifically  identified  || . 

LF.  The  power  in  the  low-frequency  range:  0.04-0.15  cycles/interval.  The  LF  band  may 
reflect  both  sympathetic  and  vagal  activity  but  its  interpretation  is  controversial  || . 


HF.  The  power  in  the  high-frequency  range:  0.15-0.4  cycles/interval.  Efferent  vagal  activity 


is  a  major  contributor  to  the  HF  band  |26,  27,  28 


LF/HF.  The  ratio  of  the  low-frequency-range  power  to  that  in  the  high-frequency  range. 
This  ratio  may  mirror  either  sympatho- vagal  balance  or  reflect  sympathetic  modulations  [^J . 
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2.3  Standard  Time-Domain  Measures 


We  consider  three  time-domain  measures  commonly  used  in  HRV  analysis.  The  first  and 
last  are  highly  correlated  with  each  other  inasmuch  as  they  estimate  the  high-frequency 
variations  in  the  heart  rate  p.  They  are: 

pNN50.  The  relative  proportion  of  successive  NN  intervals  (normal-to-normal  intervals, 
i.e. ,  all  intervals  between  adjacent  QRS  complexes  resulting  from  sinus  node  depolarizations 
IQ)  with  interval  differences  greater  than  50  ms. 

SDANN.  The  Standard  Deviation  of  the  Average  NN  interval  calculated  in  five-minute 
segments.  It  is  often  calculated  over  a  24-hour  period.  This  measure  estimates  fluctuations 
over  frequencies  smaller  than  0.003  cycles/sec. 

SDNN  ((Tint).  The  Standard  Deviation  of  the  NN  interval  set  {r*}  specified  in  units 
of  seconds.  This  measure  is  one  of  the  more  venerable  among  the  many  scale-dependent 
measures  that  have  long  been  used  for  HRV  analysis  0  §  0  0- 


2.4  Other  Standard  Measures 


There  are  several  other  well-known  measures  that  have  been  considered  for  HRV  analysis.  For 
completeness,  we  briefly  mention  two  of  them  here:  the  event-number  histogram  and  the  Fano 
factor  P  P .  Just  as  the  interevent-interval  histogram  provides  an  estimate  of  the  probability 
density  function  of  interevent-interval  magnitude,  the  event-number  histogram  provides  an 
estimate  of  the  probability  mass  function  of  the  number  of  events.  Construction  of  the  event- 
number  histogram,  like  the  interevent-interval  histogram,  involves  loss  of  information,  in  this 
case  the  ordering  of  the  counts.  However,  whereas  the  time  scale  of  information  contained 
in  the  interevent-interval  histogram  is  the  mean  interevent  interval,  which  is  intrinsic  to  the 
process  under  consideration,  the  event-number  histogram  reflects  behavior  occurring  on  the 
adjustable  time  scale  of  the  counting  window  T .  The  Fano  factor,  which  is  the  variance  of 
the  number  of  events  in  a  specified  counting  time  T  divided  by  the  mean  number  of  events 
in  that  counting  time,  is  a  measure  of  correlation  over  different  time  scales  T .  This  measure 
is  sometimes  called  the  index  of  dispersion  of  counts  [RTI] .  In  terms  of  the  sequence  of  counts 


illustrated  in  Fig.  1,  the  Fano  factor  is  simply  the  variance  of  {IVj}  divided  by  the  mean  of 

TO- 
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2.5  Novel  Scale-Dependent  Measures 


The  previous  standard  measures  are  all  well-established  scale-dependent  measures.  We  now 
describe  a  set  of  recently  devised  scale-dependent  measures  whose  performance  we  evaluate. 
Throughout  this  chapter,  when  referring  to  intervals,  we  denote  the  fixed  scale  as  m;  when 
referring  to  time,  we  employ  T. 


2.5.1  Allan  Factor  [A(T)] 


In  this  section  we  present  a  measure  we  first  defined  in  1996  [[S|  and  called  the  Allan  factor. 
We  quickly  found  that  this  quantity  was  a  useful  measure  of  HRV  [|j.  The  Allan  factor  is 
the  ratio  of  the  event-number  Allan  variance  to  twice  the  mean: 


E{[iVm(r)-iv.(r)]2} 

1  '  2E{Ni+1(T)} 


(9) 


The  Allan  variance,  as  opposed  to  the  ordinary  variance,  is  defined  in  terms  of  the  variability 
of  successive  counts  [0  S  0-  As  such,  it  is  a  measure  based  on  the  Haar  wavelet.  The 
Allan  variance  was  first  introduced  in  connection  with  the  stability  of  atomic-based  clocks 
B-  Because  the  Allan  factor  functions  as  a  derivative,  it  has  the  salutary  effect  of  mitigating 


against  linear  nonstationarities. 


The  Allan  factor  of  a  point  process  generally  varies  as  a  function  of  the  counting  time 
T;  the  exception  is  the  homogeneous  Poisson  point  process.  For  a  homogeneous  Poisson 
point  process,  A(T )  =  1  for  any  counting  time  T.  Any  deviation  from  unity  in  the  value 
of  A(T )  therefore  indicates  that  the  point  process  in  question  is  not  Poisson  in  nature.  An 
excess  above  unity  reveals  that  a  sequence  is  less  ordered  than  a  homogeneous  Poisson  point 
process,  while  values  below  unity  signify  sequences  which  are  more  ordered.  For  a  point 
process  without  overlapping  events  the  Allan  factor  approaches  unity  as  T  approaches  zero. 


A  more  complex  wavelet  Allan  factor  can  be  constructed  to  eliminate  polynomial  trends 
35],  |36],  [J7J.  The  Allan  variance,  E[(W+i  —  A*)2]  may  be  recast  as  the  variance  of  the  integral 
of  the  point  process  under  study  multiplied  by  the  following  function: 


”0Haar(t) 


—  1  for  —  T  <  t  <  0, 
+1  for  0  <  t  <  T, 

0  otherwise. 


(10) 


Equation  (]TD[)  defines  a  scaled  wavelet  function,  specifically  the  Haar  wavelet.  This  can  be 
generalized  to  any  admissible  wavelet  -0( t );  when  suitably  normalized  the  result  is  a  wavelet 
Allan  factor  IE5L 


13 


2.5.2  Wavelet- Transform  Standard  Deviation  [crwav(m)] 


Wavelet  analysis  has  proved  to  be  a  useful  technique  for  analyzing  signals  at  multiple 
scales  |h|,  fit],  (IT],  f|2|,  [|3|,  ft 


ff3[,  [h^|.  It  permits  the  time  and  frequency  characteristics  of  a  sig¬ 
nal  to  be  simultaneously  examined,  and  has  the  advantage  of  naturally  removing  polynomial 
nonstationarities  |36],  57],  ffH] .  The  Allan  factor  served  in  this  capacity  for  the  counting  process 


{ Ni },  as  discussed  above.  Wavelets  similarly  find  use  in  the  analysis  of  RR-interval  series. 
They  are  attractive  because  they  mitigate  against  the  nonstationarities  and  slow  variations 
inherent  in  the  interbeat-interval  sequence.  These  arise,  in  part,  from  the  changing  activity 
level  of  the  subject  during  the  course  of  a  24-hour  period. 


Wavelet  analysis  simultaneously  gives  rise  to  both  scale-dependent  and  scale-independent 
measures  J46|,  affording  the  experimenter  an  opportunity  to  compare  the  two  approaches. 
In  this  latter  capacity  wavelet  analysis  provides  an  estimate  of  the  wavelet-transform  fractal 
(scaling)  exponent  aw 


m 


as  discussed  in  the  context  of  HRV  in  subsections  12.6.2 


and  [3.6. 1|.  As  a  result  of  these  salutary  properties  we  devote  particular  attention  to  the 
wavelet  analysis  of  HRV  in  this  Chapter. 


A  dyadic  discrete  wavelet  transform  for  the  RR-interval  sequence  {rj}  may  be  defined  as 


inn 


m 


L- 1 


Wm,n(m)  =  —=  Y  Ti7p(i/m  -  n). 


(11) 


i=0 


The  quantity  -0  is  the  wavelet  basis  function,  and  L  is  the  number  of  RR  intervals  in  the 
set  {Tj}.  The  scale  m  is  related  to  the  scale  index  j  by  m  =  2J.  Both  j  and  the  translation 
variable  n  are  nonnegative  integers.  The  term  dyadic  refers  to  the  use  of  scales  that  are 
integer  powers  of  2.  This  is  an  arbitrary  choice;  the  wavelet  transform  could  be  calculated  at 
arbitrary  scale  values,  although  the  dyadic  scale  enjoys  a  number  of  convenient  mathematical 
properties  HTL 


m 


The  dyadic  discrete  wavelet  transform  calculated  according  to  this  prescription  generates 
a  three-dimensional  space  from  a  two-dimensional  signal  graph.  One  axis  is  time  or,  in 
our  case,  the  RR-interval  number  the  second  axis  is  the  scale  m;  and  the  third  axis  is 
the  strength  of  the  wavelet  component.  Pictorially  speaking,  the  transform  gives  rise  to  a 
landscape  whose  longitude  and  latitude  are  RR-interval  number  and  scale  of  observation, 
while  the  altitude  is  the  value  of  the  discrete  wavelet  transform  at  the  interval  i  and  the 
scale  m. 


Figure  2  provides  an  example  of  such  a  wavelet  transform,  where  ip{x)  is  the  simple  Haar 
wavelet.  Figure  2(a)  illustrates  the  original  wavelet,  a  function  that  is  by  definition  =  1 
for  x  between  0  and  0.5;  -0(x)  =  —1  for  x  between  0.5  and  1;  and  0(x)  =  0  elsewhere.  Figure 
2(b)  illustrates  the  wavelet  scaled  by  the  factor  m  =  16,  which  causes  it  to  last  for  16  samples 
rather  than  1;  and  delayed  by  a  factor  of  n  =  3  times  the  length  of  the  wavelet,  so  that  it 
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begins  at  x  —  nm  =  48.  Figure  2(c)  shows  a  sequence  of  interbeat-interval  values  multiplied 
by  the  scaled  and  shifted  wavelet  [the  summand  in  Eq.  ([□])].  The  abscissa  is  labeled  i  rather 
than  x  to  indicate  that  we  have  a  discrete-time  process  comprised  of  the  sequence  {t/}.  In 
this  particular  example,  only  values  of  t,  between  i  =  48  and  63  survive.  Adding  them  (with 
the  appropriate  sign)  provides  the  wavelet  transform  beginning  at  interval  number  i  =  48  at 
a  scale  of  m  =  16. 


For  the  Haar  wavelet  the  calculation  of  the  wavelet  transform  is  therefore  tantamount  to 
adding  the  eight  RR  intervals  between  intervals  48  and  55  inclusive,  and  then  subtracting 
the  eight  subsequent  RR  intervals  between  intervals  56  and  63  inclusive,  as  illustrated  in 
Fig.  2(c).  Moving  this  window  across  interval  number  allows  us  to  see  how  the  wavelet 
transform  evolves  with  interval  number,  whereas  varying  the  scale  of  the  window  permits 
this  variation  to  be  observed  over  a  range  of  resolutions,  from  fine  to  coarse  (smaller  scales 
allow  the  observation  of  more  rapid  variations,  i.e.  higher  frequencies). 


A  simple  measure  that  can  be  devised  from  the  wavelet  transformation  is  the  standard 


deviation  of  the  wavelet  transform  as  a  function  of  scale  46,  47,  48 


cq. 


ym)  = 


E 


{\Wm,n(rn)  -E[Wm,n(m)]\2} 


1/2 


(12) 


where  the  expectation  is  taken  over  the  process  of  RR  intervals,  and  is  independent  of  n.  It 
is  readily  shown  that  E [IFm,n(^)]  =  0  for  all  values  of  m  so  that  a  simplified  form  for  the 
wavelet-transform  standard  deviation  emerges: 


crv 


(m)  —  |e  |Wm,n(m)|2]} 


1/2 


(13) 


This  quantity  has  recently  been  shown  to  be  quite  valuable  for  HRV  analysis 
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50j.  The  special  case  obtained  by  using  the  Haar- wavelet  basis  and  evaluating  Eq.  (fL3])  at 
m  —  l  yields  the  standard  deviation  of  the  difference  between  pairs  of  consecutive  interbeat 
intervals.  This  special  case  is  therefore  identical  to  the  well-known  HRV  measure  referred  to 
as  RMSSD  [§,  an  abbreviation  for  Root-Mean-Square  of  Successive-interval  Differences. 


Figure  3  provides  an  example  in  which  the  discrete  wavelet  transform  is  calculated  using 
an  RR-interval  data  set.  In  Fig.  3(a),  the  original  RR  interbeat-interval  series  is  shown, 
while  Fig.  3(b)  shows  the  dyadic  discrete  wavelet  transform  at  three  different  scales  as  a 
function  of  RR-interval  number.  It  is  important  and  interesting  to  note  that  the  trends  and 
baseline  variations  present  in  the  original  time  series  have  been  removed  by  the  transform. 
As  illustrated  in  Fig.  3(b),  the  wavelet-transform  standard  deviation  <rwav  typically  increases 
with  the  scale  m.  When  plotted  versus  scale,  this  quantity  provides  information  about  the 
behavior  of  the  signal  at  all  scales.  In  Sec.  [|  we  show  how  this  measure  can  be  effectively 
used  to  separate  heart-failure  patients  from  normal  normal  subjects. 
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2.5.3  Relationship  of  Wavelet  [<xwav(ra)]  and  Spectral  Measures  [S  (/)] 


Is  there  a  spectral  measure  equivalent  to  the  wavelet-transform  standard  deviation?  We 
proceed  to  show  that  the  wavelet-transform  standard  deviation  <jwav  (m)  and  the  interval- 
based  power  spectral  density  ST(f)  are  isomorphic  [|49|,  so  that  the  answer  is  yes  under 
conditions  of  stationarity.  Though  their  equivalence  is  most  easily  analyzed  in  the  continuous 
domain,  the  results  are  readily  translated  to  the  discrete  domain  by  interpreting  the  discrete 
wavelet  transform  as  a  discretized  version  of  a  continuous  wavelet  transform. 


The  continuous  wavelet  transform  of  a  signal  r(t)  is  defined  as 

dt 


1  /  f  —  V 

WT(s,r)  =  —  r(t)ip* 


S  J— oo 


(14) 


where  s  and  r  are  continuous-valued  scale  and  translation  parameters  respectively,  if  is  a 
wavelet  basis  function,  and  *  denotes  complex  conjugation.  Since  E[Wr]  =  0,  the  variance 
of  WT  at  scale  s  is 

B(s)  =  <4,v(s)  =  E[|W'r(s,r)|2],  (15) 

which  can  be  written  explicitly  as 


D(s)  =  E 


S  J-c 


r{tW 


t  —  r 


dt 


s  J-c 


r*(f>  (~ - -W 


For  a  wide-sense  stationary  signal  the  variance  can  be  written  as 


DM  = 


|  roc  roc 


S  J —oo  J  —oo 


R(t  -  t'w 


t  —  r 


'  t'  —  r ' 


dtdt' 


(16) 


(V) 


where  R  is  the  autocorrelation  function  of  r.  Routine  algebraic  manipulation  then  leads  to 

/OO 

R{sy)W^,{\,y)dy  (18) 

-OO 

or,  alternatively, 


D(s)  =  s 


f  =  oo 


StU) 


iy=-o o 


1,  y)  exp(j2Tcfsy)dy 


df 


(19) 


where  W^(l,y)  is  the  wavelet  transform  of  the  wavelet  itself  (termed  the  wavelet  kernel), 
and  ST(f)  is  the  power  spectral  density  of  the  signal. 

For  the  dyadic  discrete  wavelet  transform  that  we  have  used,  Eq.  (|19|)  becomes 


roo  r  poo 

D(m)  =  oh-LvM  =  m  ST(f )  /  W^{l,y)exp(J2nfmy)dy 

J  f=— oo  LJy=—oo 

/oo 

Sr{f)H(mf)df. 

-OO 


df 


(20) 
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We  conclude  that  for  stationary  signals  the  interval-based  power  spectral  density  ST(f)  is 
directly  related  to  the  wavelet-transform  standard  deviation  crwav(m)  through  an  integral 
transform.  This  important  result  has  a  simple  interpretation:  the  factor  in  square  brackets 
in  Eq.  (^)  represents  a  bandpass  filter  H(mf )  that  only  passes  spectral  components  in  a 
bandwidth  surrounding  the  frequency  fm  that  corresponds  to  the  scale  m.  This  is  because  the 
Fourier  transform  of  a  wavelet  kernel  is  constrained  to  be  bandpass  in  nature.  For  a  discrete- 
index  sequence,  the  sampling  “time”  can  be  arbitrarily  set  to  unity  so  that  a  frequency  fm 
corresponds  to  1/m.  We  conclude  that  information  obtained  from  a  D(m)~ based  statistic  is 
also  accessible  through  interval-based  power  spectral  density  measures.  In  Sec.  ||  we  explicitly 
show  that  the  two  measures  are  comparable  in  their  abilities  to  discriminate  heart-failure 
patients  from  normal  subjects. 


2.5.4  Detrended  Fluctuation  Analysis  [DFA(m)] 


Detrended  fluctuation  analysis  (DFA)  was  originally  proposed  as  a  technique  for  quantifying 
the  nature  of  long-range  correlations  in  a  time  series  [J5T],  [)2],  Q.  As  implied  by  its  name, 
it  was  conceived  as  a  method  for  detrending  variability  in  a  sequence  of  events.  The  DFA 
computation  involves  the  calculation  of  the  summed  series 

k 

y(k)  =  I]{a*e[t]}  (21) 

Z— 1 


where  y(k )  is  the  hth  value  of  the  summed  series  and  E[r]  denotes  the  average  over  the  set 
{Vi}.  The  summed  series  is  then  divided  into  segments  of  length  m  and  a  least-squares  fit  is 
performed  on  each  of  the  data  segments,  providing  the  trends  for  the  individual  segments. 
Detrending  is  carried  out  by  subtracting  the  local  trend  ym(k)  in  each  segment.  The  root- 
mean-square  fluctuation  of  the  resulting  series  is  then 


F(m) 


k= 1 


1/2 


Vm{k )] 


(22) 


The  functional  dependence  of  F{m)  is  obtained  by  evaluations  over  all  segment  sizes  m. 


Although  detrended  fluctuation  analysis  was  originally  proposed  as  a  method  for  esti¬ 
mating  the  scale-independent  fractal  exponent  of  a  time  series  [^T],  as  discussed  in  Sec.  [2.6. ![, 
we  consider  its  merits  as  a  scale- dependent  measure.  As  will  be  demonstrated  in  Sec.  [3].  a 
plot  of  F(m)  versus  m  reveals  a  window  of  separation  between  congestive-heart-failure  pa¬ 
tients  and  normal  subjects  over  a  limited  range  of  scales,  much  as  that  provided  by  the  other 
scale-dependent  measures  discussed  in  this  section.  Because  DFA  is  an  ad  hoc  measure  that 
involves  nonlinear  computations  it  is  difficult  to  relate  it  to  other  scale-dependent  measures 
in  the  spirit  of  Eq.  (|0[).  Furthermore,  as  will  become  clear  in  Sec.  |3.6.2|,  relative  to  other 
measures  DFA  is  highly  time  intensive  from  a  computational  point-of-view. 
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2.6  Scale-Independent  Measures 


Scale-independent  measures  are  designed  to  estimate  fractal  exponents  that  characterize 
scaling  behavior  in  one  or  more  statistics  of  a  sequence  of  events,  as  discussed  in  Sec.  p.1.2 


The  canonical  example  of  a  scale-independent  measure  in  HRV  is  the  fractal  exponent  a$  of 
the  interbeat-interval  power  spectrum,  associated  with  the  decreasing  power-law  form  of  the 
spectrum  at  sufficiently  low  frequencies  f:  ST(f)  oc  /  “s  0  0  0-  Other  scale-independent 
measures  have  been  examined  by  us  |[7].  0.  |TT|.  [46|  and  by  others  |[d]  0,  |56],  [I]|  in  connection 

with  HRV  analysis.  For  exponent  values  encountered  in  HRV,  and  infinite  data  length,  all 
measures  should  in  principle  lead  to  a  unique  fractal  exponent.  In  practice,  however,  finite 
data  length  and  other  factors  introduce  bias  and  variance,  so  that  different  measures  give 
rise  to  different  results.  The  performance  of  scale-independent  measures  has  been  compared 


with  that  of  scale-dependent  measures  for  assessing  cardiac  dysfunction  [46,  47 


2.6.1  Detrended-Fluctuation- Analysis  Power-Law  Exponent  (o:D  ) 


The  DFA  technique,  and  its  use  as  a  scale-dependent  measure,  has  been  described  in 
Sec.  |2.5.4|.  A  number  of  recent  studies  [ETl  1531  |56|  have  considered  the  extraction  of  power- 


law  exponents  from  DFA  and  their  use  in  HRV.  As  originally  proposed  [[nj,  log[F(m)]  is 
plotted  against  log(m)  and  scaling  exponents  are  obtained  by  fitting  straight  lines  to  sections 
of  the  resulting  curve  -  the  exponents  are  simply  the  slopes  of  the  linearly  fitted  segments 
on  this  doubly  logarithmic  plot.  The  relationship  between  the  scaling  exponents  has  been 


proposed  as  a  means  of  differentiating  normal  from  pathological  subjects  [51 


2.6.2  Wavelet- Transform  Power-Law  Exponent  (ckw  ) 


The  use  of  the  wavelet  transform  as  a  scale-dependent  measure  was  considered  in  Sec.  |2.5.2 


It  was  pointed  out  that  a  scale-independent  measure  also  emerges  from  the  wavelet-transform 
standard  deviation.  The  wavelet-transform  fractal  exponent  aw  is  estimated  directly  from 
the  wavelet  transform  as  twice  the  slope  of  the  curve  log[awav(m)]  versus  log(m),  measured 
at  large  values  of  m  [[46| .  The  factor  of  two  is  present  because  the  fractal  exponent  is  related 
to  variance  rather  than  to  standard  deviation. 


2.6.3  Periodogram  Power-Law  Exponent  (cks) 

The  description  of  the  periodogram  as  a  scale-dependent  measure  was  provided  in  Sec.  [2~2l 
The  periodogram  fractal  exponent  as  0  0]  is  obtained  as  the  least-squares-fit  slope  of 
the  spectrum  when  plotted  on  doubly  logarithmic  coordinates.  The  range  of  low  frequencies 
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over  which  the  slope  is  estimated  stretches  between  10/L  and  100/L  where  L  is  the  length 
of  the  data  set  [|17] . 


2.6.4  Allan-Factor  Power-Law  Exponent  (o:A  ) 


The  use  of  the  Allan  factor  as  a  scale-dependent  measure  was  considered  in  Sec.  [2.5. 1|.  The 
Allan  factor  fractal  exponent  a  a  ||  0  is  obtained  by  determining  the  slope  of  the  best¬ 
fitting  straight  line,  at  large  values  of  T,  to  the  Allan  factor  curve  [Eq.  @]  plotted  on  doubly 
logarithmic  coordinates.  Estimates  of  a  obtained  from  the  Allan  factor  can  range  up  to  a 
value  of  three  [32],  The  use  of  wavelets  more  complex  than  the  Haar  enables  an  increased 


range  of  fractal  exponents  to  be  accessed,  at  the  cost  of  a  reduction  in  the  range  of  counting 
time  over  which  the  wavelet  Allan  factor  varies  as  TaA.  In  general,  for  a  particular  wavelet 
with  regularity  (number  of  vanishing  moments)  R ,  fractal  exponents  a  <  2 R  +  1  can  be 
reliably  estimated  ]36],  |38| .  For  the  Haar  basis,  R  =  1  whereas  all  other  wavelet  bases  have 
R  >  1.  A  wavelet  Allan  factor  making  use  of  bases  other  than  the  Haar  is  therefore  required 
for  fractal-rate  stochastic  point  processes  for  which  a  >  3.  For  processes  with  a  <  3, 
however,  the  Allan  factor  appears  to  be  the  best  choice  pi  RSI . 


2.6.5  Rescaled-Range- Analysis  Power-Law  Exponent  (aR  ) 


Rescaled  range  analysis  [[19],  [77].  [79]| ,  provides  information  about  correlations  among  blocks 

of  interevent  intervals.  For  a  block  of  k  interevent  intervals,  the  difference  between  each 
interval  and  the  mean  interevent  interval  is  obtained  and  successively  added  to  a  cumulative 
sum.  The  normalized  range  R(k )  is  the  difference  between  the  maximum  and  minimum 
values  that  the  cumulative  sum  attains,  divided  by  the  standard  deviation  of  the  interval 
size.  R[k)  is  plotted  against  k.  Information  about  the  nature  and  the  degree  of  correlation 
in  the  process  is  obtained  by  fitting  R[k)  to  the  function  kH ,  where  H  is  the  so-called 
Hurst  exponent  0-  For  H  >  0.5  positive  correlation  exists  among  the  intervals,  whereas 
H  <  0.5  indicates  the  presence  of  negative  correlation;  H  =  0.5  obtains  for  intervals  with 
no  correlation.  Renewal  processes  yield  H  =  0.5.  For  negatively  correlated  intervals,  an 
interval  that  is  larger  than  the  mean  tends,  on  average,  to  be  preceded  or  followed  by  one 
smaller  than  the  mean. 


The  Hurst  exponent  H  is  generally  assumed  to  be  well  suited  to  processes  that  exhibit 
long-term  correlation  or  have  a  large  variance 


but  there  are  limits  to  its 
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58 


robustness  since  it  exhibits  large  systematic  errors  and  highly  variable  estimates  for  some 

m 


fractal  sequences  ||  |60],  [H|.  Nevertheless,  it  provides  a  useful  indication  of  correlation  in  a 
sequence  of  events  arising  from  the  ordering  of  the  interevent  intervals  alone. 

The  exponent  is  ambiguously  related  to  the  Hurst  exponent  H ,  since  some  authors 
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have  used  the  quantity  H  to  index  fractal  Gaussian  noise  whereas  others  have  used  the  same 
value  of  H  to  index  the  integral  of  fractal  Gaussian  noise  (which  is  fractional  Brownian 
motion).  The  relationship  between  the  quantities  is  =  2 H  —  1  for  fractal  Gaussian  noise 
and  oir  =  2 H  +  1  for  fractal  Brownian  motion.  In  the  context  of  this  work,  the  former 
relationship  holds. 


2.7  Estimating  the  Performance  of  a  Measure 


We  have,  to  this  point,  outlined  a  variety  of  candidate  measures  for  use  in  HRV  analysis. 
The  task  now  is  to  determine  the  relative  value  of  these  measures  from  a  clinical  perspective. 
We  achieve  this  by  turning  to  estimation  theory  [|62  . 


A  statistical  measure  obtained  from  a  finite  set  of  actual  data  is  characterized  by  an 
estimator.  The  fidelity  with  which  the  estimator  can  approximate  the  true  value  of  the 
measure  is  determined  by  its  bias  and  variance.  The  bias  is  the  deviation  of  the  expected 
value  of  the  estimator  from  its  true  underlying  value  (assuming  that  this  exists)  whereas  the 
variance  indicates  the  expected  deviation  from  the  the  mean.  An  ideal  estimator  has  zero 
bias  and  zero  variance,  but  this  is  not  achievable  with  a  finite  set  of  data.  For  any  unbiased 
estimator  the  Cramer-Rao  bound  provides  a  lower  bound  for  the  estimator  variance;  measures 
that  achieve  the  Cramer-Rao  bound  are  called  efficient  estimators.  The  estimator  bias  and 
variance  play  a  role  in  establishing  the  overall  statistical  significance  of  conclusions  based  on 
the  value  returned  by  the  estimator. 


2.7.1  Statistical  Significance:  p,  d°,  h,  and  d 


The  concept  of  statistical  significance  extends  the  basic  properties  of  bias  and  variance  [33 


It  provides  a  probabilistic  interpretation  of  how  likely  it  is  that  a  particular  value  of  the 
estimator  might  occur  by  chance  alone,  arising  from  both  random  fluctuations  in  the  data 
and  the  inherent  properties  of  the  estimator  itself. 


A  frequently  used  standard  of  statistical  significance  is  the  p-value,  the  calculation  of 
which  almost  always  implicitly  assumes  a  Gaussian-distributed  dependent  variable.  A  lower 
p-value  indicates  greater  statistical  significance,  and  a  measure  is  said  to  be  statistically 
significant  to  a  value  of  po  when  p  <  p$.  The  distributions  obtained  from  HRV  measurements 
are  generally  not  Gaussian,  however,  so  that  the  usual  method  for  estimating  the  p-value 
cannot  be  used  with  confidence.  Since  other  methods  for  estimating  the  p- value  require  more 
data  than  is  available  we  do  not  consider  this  quantity  further. 

Another  often-used  distribution-dependent  standard  is  the  d'-value.  It  serves  to  indicate 
the  degree  of  separation  between  two  distributions,  and  has  been  widely  used  in  signal 
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detection  theory  and  psychophysics  where  the  two  distributions  represent  noise  and  signal- 
plus-noise  [Q.  The  most  common  definition  of  df  is  the  difference  in  the  means  of  two 
Gaussian  distributions  divided  by  their  common  standard  deviation.  Two  closely  related 
distribution-dependent  cousins  of  d'  are  the  detection  distance  h,  defined  as  the  difference 
in  the  means  of  the  two  Gaussian  distributions  divided  by  the  square-root  of  the  sum  of 
their  variances;  and  the  detection  distance  d,  defined  as  the  difference  in  the  means  of  the 
two  Gaussian  distributions  divided  by  the  sum  of  their  standard  deviations.  Larger  values 
of  d1,  h,  and  d  indicate  improved  separation  between  the  two  distributions  and  therefore 
reduced  error  in  assigning  an  outcome  to  one  or  the  other  of  the  hypotheses. 


Because  HRV  measures  are  intended  to  provide  diagnostic  information  in  a  clinical 
setting,  and  do  not  return  Gaussian  statistics,  the  evaluation  of  their  performance  using 
distribution-independent  means  is  preferred.  Two  techniques  for  achieving  this,  positive  and 
negative  predictive  values,  and  receiver  operator  characteristic  (ROC)  analysis,  are  described 
below.  Neither  requires  knowledge  of  the  statistical  distribution  of  the  measured  quantities 
and  both  are  useful. 


2.7.2  Positive  and  Negative  Predictive  Values 


The  performance  of  the  various  HRV  measures  discussed  previously  can  be  effectively  com¬ 
pared  using  positive  predictive  values  and  negative  predictive  values,  the  proportion  of  cor¬ 
rect  positive  and  negative  identifications  respectively.  When  there  is  no  false  positive  (or 
negative)  detection,  the  predictive  value  is  equal  to  unity  and  there  is  perfect  assignment. 
Furthermore,  when  the  individual  values  of  a  measure  for  normal  subjects  and  patients  do 
not  overlap,  the  predictive  value  curves  are  typically  monotonic,  either  increasing  or  decreas¬ 
ing,  with  the  threshold.  A  detailed  discussion  of  positive  and  negative  predictive  values  is 


provided  in  Sec.  3.4 


2.7.3  Receiver-Operating-Characteristic  (ROC)  Analysis 


Receiver-operating-characteristic  (ROC)  analysis  is  an  objective  and  highly  effec¬ 

tive  technique  for  assessing  the  performance  of  a  measure  when  it  is  used  in  binary  hypothesis 
testing.  This  format  provides  that  a  data  sample  be  assigned  to  one  of  two  hypotheses  or 
classes  (e.g.,  pathologic  or  normal)  depending  on  the  value  of  some  measured  statistic  relative 
to  a  threshold  value.  The  efficacy  of  a  measure  is  then  judged  on  the  basis  of  its  sensitivity 
(the  proportion  of  pathologic  patients  correctly  identified)  and  its  specificity  (the  propor¬ 
tion  of  normal  subjects  correctly  identified).  The  ROC  curve  is  a  graphical  presentation  of 
sensitivity  versus  1— specificity  as  a  threshold  parameter  is  swept.  Note  that  sensitivity  and 
specificity  relate  to  the  status  of  the  patients  (pathologic  and  normal)  whereas  predictive 
values  relate  to  the  status  of  the  identifications  (positive  and  negative). 
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The  area  under  the  ROC  curve  serves  as  a  well-established  index  of  diagnostic  accuracy 
|G4|.  |65|;  the  maximum  value  of  1.0  corresponds  to  perfect  assignment  (unity  sensitivity  for 
all  values  of  specificity)  whereas  a  value  of  0.5  arises  from  assignment  to  a  class  by  pure 
chance  (areas  <  0.5  arise  when  the  sense  of  comparison  is  reversed).  ROC  analysis  can  be 
used  to  choose  the  best  of  a  host  of  different  candidate  diagnostic  measures  by  comparing 
their  ROC  areas,  or  to  establish  for  a  single  measure  the  tradeoff  between  data  length  and 
misidentihcations  (misses  and  false  positives)  by  examining  ROC  area  as  a  function  of  record 
length.  A  minimum  record  length  can  then  be  specified  to  achieve  acceptable  classification 
accuracy. 


As  pointed  out  above,  ROC  analysis  relies  on  no  implicit  assumptions  about  the  statistical 
nature  of  the  data  set  |65|,  so  that  it  is  generally  more  suitable  [[47]  for  analyzing  non- 
Gaussian  time  series  than  are  measures  of  statistical  significance  such  as  p- value,  h,  and  d. 
Another  important  feature  of  ROC  curves  is  that  they  are  insensitive  to  the  units  employed 
(e.g.,  spectral  magnitude,  magnitude  squared,  or  log  magnitude);  ROC  curves  for  a  measure 
M  are  identical  to  those  for  any  monotonic  transformation  thereof  such  as  Mx  or  log(M). 
In  contrast  the  values  of  d',  h,  and  d  are  generally  modified  by  such  transformations,  as  will 


be  demonstrated  in  Sec.  13.5.1. 


3  Discriminating  Heart-Failure  Patients  from  Normal 
Subjects 


We  now  proceed  to  examine  the  relative  merits  of  various  HRV  measures  for  discriminating 
congestive- heart-failure  (CHF)  patients  from  normal  subjects.  Specifically  we  contrast  and 
compare  the  performance  of  the  16  measures  set  forth  in  Sec.  [2|  VLF,  LF,  HF,  LF/HF, 
pNN50,  SDANN,  SDNN  (<jint),  A(T ),  crwav(m),  ST(f ),  DFA(m),  an,  ctw,  T4,  and  a.R. 

After  discussing  the  selection  of  an  appropriate  scale  m,  we  use  predictive  value  plots 
and  ROC  curves  to  select  a  particular  subset  of  HRV  markers  that  appears  to  be  promising 
for  discerning  the  presence  of  heart  failure  in  a  patient  population. 


3.1  Database 

The  RR  recordings  analyzed  in  this  section  were  drawn  from  the  Beth-Israel  Hospital  (Boston, 
MA)  Heart-Failure  Database  which  includes  12  records  from  normal  subjects  (age  29-64 
years,  mean  44  years)  and  12  records  from  severe  congestive-heart-failure  patients  (age  22- 
71  years,  mean  56  years).  The  recordings  were  made  with  a  Holter  monitor  digitized  at  a 
fixed  value  of  250  samples/sec.  Also  included  in  this  database  are  3  RR  records  for  CHF  pa- 
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tients  who  also  suffered  from  atrial  fibrillation  (AF);  these  records  are  analyzed  as  a  separate 
class.  All  records  contain  both  diurnal  and  nocturnal  segments.  The  data  were  originally 
provided  to  us  in  1992  by  D.  Rigney  and  A.  L.  Goldberger. 

A  detailed  characterization  of  each  of  the  records  is  presented  in  Table  1  of  Ref.  [|| ;  some 
statistical  details  are  provided  in  Table  Al.  Of  the  27  recordings,  the  shortest  contained 
Lm™  =  75821  RR  intervals;  the  remaining  26  recordings  were  truncated  to  this  length  before 
calculating  the  16  HRV  measures. 


3.2  Selecting  a  Scale 


A  value  for  the  scale  m  that  suitably  discriminates  heart-failure  patients  from  normal  subjects 
can  be  inferred  from  our  recent  wavelet  studies  of  the  CHF  and  normal  records  from  the 
same  database  as  discussed  in  Sec.  |3.1|  [47|.  With  the  help  of  the  wavelet-transform 


standard  deviation  crwav(m)  discussed  in  detail  in  Sec.  [2.5.2|,  we  discovered  a  critical  scale 
window  near  m  =  32  interbeat  intervals  over  which  the  normal  subjects  exhibited  greater 
fluctuations  than  those  afflicted  with  heart  failure.  For  these  particular  long  data  sets,  we 
found  that  it  was  possible  to  perfectly  discriminate  between  the  two  groups  [|46l  |7l  |8l . 


The  results  are  displayed  in  Fig.  4,  where  crwav(m)  is  plotted  vs.  wavelet  scale  m  for  the 
12  normal  subjects  (+),  the  12  CHF  patients  s  atrial  fibrillation  (x),  and  the  3  CHF  patients 
c  atrial  fibrillation  (A),  using  Haar- wavelet  analysis.  The  AF  patients  (A)  typically  fell  near 
the  high  end  of  the  non-AF  patients  (x),  indicating  greater  RR  fluctuations,  particularly  at 
small  scales.  This  results  from  the  presence  of  non-sinus  beats.  Nevertheless  it  is  evident 
from  Fig.  4  that  the  wavelet  measure  crwav  serves  to  completely  separate  the  normal  subjects 
from  the  heart-failure  patients  (both  s  and  c  AF)  at  scales  of  16  and  32  heartbeat  intervals, 
as  reported  in  Ref.  [fhl .  One  can  do  no  better.  This  conclusion  persists  for  a  broad  range 


of  analyzing  wavelets,  from  Daubechics  2-tap  (Haar)  to  Daubechies  20-tap  [R6 


The  importance  of  this  scale  window  has  been  recently  confirmed  in  an  Israeli-Danish 
study  of  diabetic  patients  who  had  not  yet  developed  clinical  signs  of  cardiovascular  disease 
H.  The  reduction  in  the  value  of  the  wavelet-transform  standard  deviation  erwav(32)  that 


leads  to  the  scale  window  occurs  not  only  for  CHF  (s  and  c  AF)  and  diabetic  patients, 
but  also  for  heart-transplant  patients  [fl8|,  [5(J ,  and  also  in  records  preceding  sudden  cardiac 
death  |46],  f||].  The  depression  of  (jwav(32)  at  these  scales  is  likely  associated  with  the  im¬ 


pairment  of  autonomic  nervous  system  function.  Baroreflex  modulations  of  the  sympathetic 
or  parasympathetic  tone  typically  lie  in  the  range  0.04-0.09  cycles/sec  (11-25  sec),  which 
corresponds  to  the  scale  where  crwav(m)  is  reduced. 

These  studies,  in  conjunction  with  our  earlier  investigations  which  revealed  a  similar 
critical  scale  window  in  the  counting  statistics  of  the  heartbeat  [R  B5J  (as  opposed  to  the 
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time- interval  statistics  under  discussion),  lead  to  the  recognition  that  scales  in  the  vicinity 
of  m  =  32  enjoy  a  special  status.  Those  measures  that  depend  on  a  particular  scale  are 
therefore  evaluated  at  m  =  32  and  /  =  1/32  in  the  expectation  that  these  values  maximize 
discriminability  in  the  more  usual  situation  when  the  two  classes  of  data  cannot  be  fully 
separated. 


3.3  Individual  Value  Plots 

Having  devised  a  suitable  scale  value  m  we  now  proceed  to  evaluate  the  16  measures  for 
all  27  normal  and  CHF  data  sets,  each  comprising  75821  RR  intervals.  The  results  are 
presented  in  Fig.  5  where  each  of  the  16  panels  represents  a  different  measure.  For  each 
measure  the  individual  values  for  normal  subjects  (+),  CHF  patients  s  AF  (x),  and  CHF 
patients  c  AF  (A)  comprise  the  left  three  columns,  respectively.  Values  in  the  right  four 
columns  correspond  to  other  cardiovascular  pathologies  and  will  be  discussed  in  Sec.  [E| 

To  illustrate  how  particular  measures  succeed  (or  fail  to  succeed)  in  distinguishing  be¬ 
tween  CHF  patients  and  normal  subjects,  we  focus  in  detail  on  two  measures:  VLF  power 
and  pNN50.  For  this  particular  collection  of  patients  and  record  lengths,  the  normal  sub¬ 
jects  all  exhibit  larger  values  of  VLF  power  than  do  the  CHF  patients;  indeed  a  horizontal 
line  drawn  at  VLF=  0.000600  completely  separates  the  two  classes.  On  the  other  hand  for 
pNN50,  though  the  normals  still  have  larger  values  on  average,  there  is  a  region  of  overlap 
of  CHF  patients  and  normal  subjects  near  0.05,  indicating  that  the  two  classes  of  patients 
cannot  be  entirely  separated  using  this  measure.  Thus  for  the  full  data  set,  comprising  75821 
RR  intervals,  VLF  succeeds  in  completely  distinguishing  CHF  patients  and  normal  subjects 
whereas  pNN50  does  not. 

Examining  all  16  panels,  we  find  that  six  measures  manage  to  completely  separate  the 
normal  subjects  (first  column)  from  the  heart-failure  patients  (second  and  third  columns) 
while  the  remaining  10  fail  to  do  so.  The  six  successful  measures  are  highlighted  by  boldface 
font  in  Fig.  5:  VLF.  LF.  A(10),  <xwav(32),  S  (1/32),  and  DFA(32). 


3.4  Predictive  Value  Plots 

How  can  the  ability  of  a  measure  to  separate  two  classes  of  subjects  be  quantified?  Returning 
to  the  VLF  panel  in  Fig.  5  we  place  a  threshold  level  9  at  an  arbitrary  position  on  the 
ordinate,  and  consider  only  the  leftmost  two  columns:  normal  subjects  and  heart-failure 
patients  who  do  not  suffer  from  atrial  fibrillation.  We  then  classify  all  subjects  for  whom 
the  VLF  values  are  <6  as  CHF  patients  (positive)  and  those  for  whom  the  VLF  values  >  6 
as  normal  (negative).  (Measures  that  yield  smaller  results  for  normal  patients,  on  average, 
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obey  a  reversed  decision  criterion.) 


If  a  subject  labeled  as  a  CHF  patient  is  indeed  so  afflicted,  then  this  situation  is  referred 
to  as  a  true  positive  (Pf)',  a  normal  subject  erroneously  labeled  as  a  CHF  patient  is  referred 
to  as  a  false  positive  (Pf)-  We  define  negative  outcomes  that  are  true  ( Nt )  and  false  (Np) 
in  an  analogous  manner.  As  pointed  out  in  Sec.  [2.7.2|,  the  positive  predictive  value  Vp  = 
Pt/(Pt  +  Pf)  and  negative  predictive  value  Vn  =  Np/(Np  +  Np)  represent  the  proportion 
of  positives  and  negatives,  respectively,  that  are  correctly  identified.  This  determination  is 
carried  out  for  many  values  of  the  threshold  9. 


Figure  6  shows  the  positive  (solid  curves)  and  negative  (dotted  curves)  predictive  values 
for  all  16  measures,  plotted  against  the  threshold  9 ,  each  in  its  own  panel.  These  curves 
are  constructed  using  the  12  normal  and  12  heart-failure  (s  AF)  records  that  comprise  the 
CHF  database  discussed  in  Sec.  ED  For  the  VLF  measure,  both  predictive  values  are 
simulataneously  unity  in  the  immediate  vicinity  of  6  =  0.000600.  This  occurs  because  Pp 
and  Np  are  both  zero  at  this  particular  value  of  9  and  reconfirms  that  the  two  classes  of 
data  separate  perfectly  in  the  VLF  panel  of  Fig.  5  at  this  threshold. 


For  threshold  values  outside  the  range  0.000544  <  9  <  0.000603,  some  of  the  patients 
will  be  incorrectly  identified  by  the  VLF  measure.  If  we  set  9  =  0.000100,  for  example,  six  of 
the  twelve  CHF  patients  will  be  incorrectly  identified  as  normal  subjects,  which  is  confirmed 
by  examining  the  VLF  panel  in  Fig.  5.  This  yields  Vn  =  Np/(Np  +  Np)  =  12/(12  +  6)  = 
0.67  <  1,  which  is  the  magnitude  of  the  negative  predictive  value  (dotted  curve)  in  the  VLF 
panel  in  Fig.  6.  At  this  value  of  the  threshold  ( 9  =  0.000100)  the  positive  predictive  value 
remains  unity  because  Pf  remains  zero. 


The  pNN50  panel  in  Fig.  5,  in  contrast,  reveals  a  range  of  overlap  in  the  individual 
values  of  the  normal  subjects  and  CHF  patients.  Consequently  as  9  increases  into  the 
overlap  region,  Vp  decreases  below  unity  and  this  happens  before  Vn  attains  unity  value. 
Thus  there  is  no  threshold  value,  or  range  of  threshold  values,  for  which  the  positive  and 
negative  predictive  values  in  Fig.  6  are  both  unity.  The  best  threshold  for  this  measure  lies 
in  the  range  0.026  <  9  <  0.050,  with  the  choice  depending  on  the  relative  benefit  of  being 
able  to  accurately  predict  the  presence  or  absence  of  CHF  in  a  patient. 

There  are  six  measures  in  Fig.  6  (indicated  in  boldface  font)  for  which  the  positive  and 
negative  predictive  values  are  both  unity  over  the  same  range  of  threshold  values.  These 
measures  are,  of  course,  the  same  six  measures  for  which  the  normal  subjects  and  heart- 
failure  patients  fall  into  disjoint  sets  in  Fig.  5. 
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3.5  ROC  Curves 


Two  other  important  clinically  relevant  quantities  that  depend  on  the  threshold  9  are  the 
sensitivity,  the  proportion  of  heart-failure  patients  that  are  properly  identified  [Pt/{Pt  + 
Np)]',  and  the  specificity,  the  proportion  of  normal  subjects  that  are  properly  identified 
[Nt/(Nt  +  Pf)]-  As  pointed  out  in  subsections  |2.7.2|  and  [2.7.3|,  sensitivity  and  specificity 
relate  to  patient  status  (pathologic  and  normal,  respectively)  whereas  the  positive  and  neg¬ 
ative  predictive  values  relate  to  identification  status  (positive  and  negative,  respectively). 
Sensitivity  and  specificity  are  both  monotonic  functions  of  the  threshold,  but  this  is  not 
generally  true  for  the  predictive  values.  The  monotonicity  property  is  salutary  in  that  it 
facilitates  the  use  of  a  parametric  plot  which  permits  these  quantities  to  be  represented  in 
compact  form.  A  plot  of  sensitivity  versus  1— specificity,  traced  out  using  various  values  of 
the  threshold  9 ,  forms  the  receiver-operating-characteristic  (ROC)  curve  (see  Sec.  p.7.3|). 


ROC  curves  are  presented  in  Fig.  7  for  all  16  measures,  again  using  the  same  12  normal 


and  12  heart-failure  (s  AF)  records  that  comprise  the  CHF  database  discussed  in  Sec.  [FT 


Because  of  the  complete  separation  between  the  two  classes  of  patients  (leftmost  two  columns 
of  the  VLF  panel  in  Fig.  5)  near  9  =  0.000600,  the  VLF  ROC  curve  in  Fig.  7  simultaneously 
achieves  unity  (100%)  sensitivity  and  unity  (100%)  specificity  (the  point  at  upper  left  corner 
of  the  ROC  curve).  For  the  pNN50  statistic,  in  contrast,  the  overlap  evident  in  Fig.  5 
prevents  this,  so  that  the  upper  left  corner  of  the  pNN50  ROC  curve  in  Fig.  7  instead  reveals 
smaller  simultaneous  values  of  sensitivity  and  specificity. 


Six  measures  in  Fig.  7  simultaneously  exhibit  unity  sensitivity  and  specificity;  these  are 
indicated  by  boldface  font  and  have  ROC  curves  that  are  perfectly  square.  They  are  clearly 
the  same  measures  for  which  the  normal  subjects  and  heart-failure  patients  fall  into  disjoint 
sets  in  Fig.  5,  and  for  which  simultaneous  positive  and  negative  predictive  values  of  unity 
are  observed  in  Fig.  6. 


3.5.1  Comparison  with  Detection-Distance  Measures 

For  didactic  purposes  we  compare  the  ROC  results  presented  immediately  above  with  those 
obtained  using  detection-distance  analysis.  As  we  indicated  in  Sec.  |2.7.1|,  care  must  be  exer¬ 
cised  when  using  these  techniques  for  anything  other  than  Gaussian-distributed  quantities. 
The  calculations  were  carried  out  using  the  same  12  normal-subject  and  12  CHF-patient 
records,  each  comprising  Lmax  =  75821  intervals.  In  Table  1  we  provide  the  detection  dis¬ 
tances  h  and  d,  in  order  of  descending  value  of  h,  for  all  16  measures.  Large  values  are  best 
since  they  indicate  that  the  two  distributions  are  well  separated. 

Five  of  the  six  measures  that  entirely  separate  the  CHF  patients  and  normal  subjects 
using  ROC  analysis  fall  in  the  top  five  positions  in  Table  1.  The  sixth  measure,  LF,  falls  in 
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ninth  position.  This  confirms  that  detection-distance  analysis  applied  to  these  long  record¬ 
ings  provides  results  that  qualitatively  agree  with  those  obtained  using  ROC  analysis.  How¬ 
ever,  detection-distance  analysis  does  not  provide  any  indication  of  how  many  (or  indeed 
whether  any)  of  the  measures  at  the  top  of  the  list  completely  separate  the  two  classes  of 
patients,  nor  does  it  provide  estimates  of  sensitivity  and  specificity.  Moreover,  the  rankings 
according  to  d  differ  from  those  according  to  h. 

Finally,  the  detection-distance  for  a  particular  measure,  as  well  as  the  relative  ranking 
of  that  measure,  depends  on  what  appear  to  be  insignificant  details  about  the  specific  form 
in  which  the  measure  is  cast.  For  example  the  h  values  for  awav(32)  and  its  square  <r^av(32) 
are  substantially  different,  and  so  are  their  rankings  in  Table  1.  As  discussed  in  Sec.  EZ3 
ROC  analysis  is  invariant  to  monotonic  transformations  of  the  measure  and  therefore  does 
not  suffer  from  this  disadvantage. 


3.5.2  ROC-Area  Curves 


Perfectly  square  ROC  curves,  associated  with  the  group  of  six  boldface-labeled  measures 
in  Figs.  5-8,  exhibit  unity  area.  These  ROCs  represent  100%  sensitivity  for  all  values  of 
specificity,  indicating  that  every  patient  is  properly  assigned  to  the  appropriate  status:  heart- 
failure  or  normal.  Though  the  perfect  separation  achieved  by  these  six  measures  endorses 
them  as  useful  diagnostic  statistics,  the  results  of  most  studies  are  seldom  so  clear-cut.  ROC 
area  will  surely  decrease  as  increasing  numbers  of  out-of-sample  records  are  added  to  the 
database,  since  increased  population  size  means  increased  variability  [|B  . 


ROC  area  also  decreases  with  diminishing  data  length;  shorter  records  yield  less  informa¬ 
tion  about  patient  condition  and  these  patients  are  therefore  more  likely  to  be  misclassified 
47l  KJ.  The  ROC  curves  in  Fig.  7  have  been  constructed  from  Holter-monitor  records  that 


contain  many  hours  of  data  (75821  RR  intervals).  It  would  be  useful  in  a  clinical  setting 
to  be  able  to  draw  inferences  from  HRV  measures  recorded  over  shorter  times,  say  minutes 
rather  than  hours.  It  is  therefore  important  to  examine  the  performance  of  the  16  HRV 
measures  as  data  length  is  decreased.  As  indicated  in  Sec.  [2.7.3|,  ROC-analysis  provides  an 


ideal  method  for  carrying  out  this  task  [47,  48 


In  Fig.  8  we  present  ROC-area  curves  as  a  function  of  the  number  of  RR  intervals  (data 
length)  L  analyzed  (64  <  L  <  75821).  The  ROC  areas  for  the  full-length  records  (Lmax  = 
75821),  which  correspond  to  the  areas  under  the  ROC  curves  presented  in  Fig.  7,  are  the 
rightmost  points  in  the  ROC-area  curves  shown  in  Fig.  8.  Results  for  shorter  records  were 
obtained  by  dividing  the  12  normal  and  12  heart-failure  (s  AF)  records  that  comprise  the 
CHF  database  (Sec.  |3. 1|)  into  smaller  segments  of  length  L.  The  area  under  the  ROC  curve 
for  that  data  length  L  was  computed  for  the  first  such  segment  for  all  16  measures,  and  then 
for  the  second  segment,  and  so  on,  for  all  segments  of  length  L  (remainders  of  length  <  L  of 
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the  original  files  were  not  used  for  the  ROC-area  calculations).  ^Frorn  the  Lmax/L  values  of 
the  ROC  area,  the  mean  and  standard  deviation  were  computed  and  plotted  in  Fig.  8.  The 
lengths  L  examined  ranged  from  L  =  26  =  64  to  L  =  216  =  65536  RR  intervals,  in  powers  of 
two,  in  addition  to  the  entire  record  of  Lmax  =  75821  intervals. 

To  illustrate  the  information  provided  by  these  curves,  we  direct  our  attention  to  the 
VLF  and  pNN50  panels  in  Fig.  8.  For  the  full-length  records  the  right-most  point  in  the 
VLF  panel  reveals  unity  area  while  that  for  pNN50  lies  somewhat  lower,  as  expected  from 
the  corresponding  ROC  curves  in  Fig.  7.  VLF  clearly  outperforms  pNN50.  As  the  data 
length  analyzed  decreases,  so  too  do  the  ROC  areas  for  both  measures  while  their  variances 
increase,  also  as  expected.  However,  when  the  data  length  dips  to  256  or  fewer  RR  intervals, 
the  performance  of  the  two  measures  reverses  so  that  pNN50  outperforms  VLF.  There  is  an 
important  point  to  be  drawn  from  this  example.  Not  only  does  the  performance  of  a  measure 
depend  on  data  length,  but  so  too  does  the  relative  performance  of  different  measures. 


3.6  Comparing  the  Measures  for  Various  Data  Lengths 

Based  on  their  overall  ability  to  distinguish  between  CHF  patients  and  normal  subjects  over 
a  range  of  data  lengths,  the  sixteen  measures  shown  in  Fig.  8  divide  roughly  into  three 
classes.  The  six  measures  that  fall  in  the  first  class,  comprising  VLF,  LF,  A(10),  <rwav(32), 
5V (1/32),  and  DFA(32),  succeed  in  completely  separating  the  two  classes  of  patients  for  data 
lengths  down  to  L  =  215  =  32768  RR  intervals.  These  six  measures  share  a  dependence  on  a 
single  scale,  or  small  range  of  scales,  near  32  heartbeat  intervals.  For  this  collection  of  data 
sets,  this  scale  appears  to  yield  the  best  performance.  Members  of  this  class  outperform  the 
other  ten  measures  at  nearly  all  data  lengths.  Apparently  the  scale  value  itself  is  far  more 
important  than  the  measure  used  to  evaluate  it. 

The  second  class,  consisting  of  HF,  the  ratio  LF/HF,  pNN50,  and  crint ,  fail  to  achieve 
complete  separation  for  any  data  size  examined.  Nevertheless  the  members  of  this  class  are 
not  devoid  of  value  in  separating  CHF  patients  from  normal  subjects.  Interestingly,  all  but 
LF/HF  provide  better  results  than  A(10),  a  member  of  the  first  class,  for  the  shortest  data 
lengths.  Results  for  these  four  measures  varied  relatively  little  with  data  size,  thus  exhibiting 
a  form  of  robustness. 

Members  of  the  third  class,  consisting  of  SDANN  and  the  five  scale-independent  measures 
«[.],  exhibit  poor  performance  at  all  data  lengths.  These  six  measures  require  long  sequences 
of  RR  intervals  to  make  available  the  long-term  fluctuations  required  for  accurate  estimation 
of  the  fractal  exponent.  Data  lengths  L  <  5000  RR  intervals  lead  to  large  variance  and 
(negative)  bias,  and  are  not  likely  to  be  meaningful.  As  an  example  of  the  kind  of  peculiarity 
that  can  emerge  when  attempting  to  apply  scale-independent  measures  to  short  records,  the 
oiA  ROC  area  decreases  below  0.5  when  the  data  size  falls  below  2048  intervals  (reversing 


the  sense  of  the  comparison  only  for  these  data  sizes  increases  the  ROC  area,  though  not 
above  0.7;  however  this  clearly  violates  the  spirit  of  the  method).  SDANN  requires  several 
5-minute  segments  to  accurately  determine  the  standard  deviation. 


3.6.1  Scale-Independent  vs  Scale-Dependent  Measures 


As  indicated  in  the  previous  subsection,  all  five  scale-independent  measures  (cud,  aw, 
a  a,  and  a#)  perform  poorly  at  all  data  lengths.  These  fractal-exponent  estimators  return 
widely  differing  results  as  is  plainly  evident  in  Fig.  5.  This  suggests  that  there  is  little  merit 
in  the  concept  of  a  single  exponent  for  characterizing  the  human  heartbeat  sequence,  no  less 
a  “universal”  one  as  some  have  proposed  |5T],  [5^,  Q. 


A  variation  on  this  theme  is  the  possibility  that  pairs  of  fractal  exponents  can  provide 
a  useful  HRV  measure.  At  small  scales  m,  Fig.  4  reveals  that  heart-failure  patients  exhibit 
smaller  values  of  the  wavelet-transform  standard-deviation  slope  than  do  normal  subjects. 
Following  Peng  et  al.  [[51 1,  who  constructed  a  measure  based  on  differences  of  DFA  scaling 
exponents  in  different  scaling  regions  in  an  attempt  to  discriminate  CHF  patients  from 
normal  subjects,  Thurner  et  al.  (46]  constructed  a  measure  based  on  differences  in  the 
wavelet-transform  standard-deviation  slope  at  different  scales.  However  the  outcome  was 
found  to  be  unsatisfactory  when  compared  with  other  available  measures;  we  concluded  the 
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same  about  the  results  obtained  by  Peng  et  al. 


Using  ROC  analysis,  as  described  in 


Sec.  [2.7.3|,  we  determined  that  the  ROC  area  for  the  measure  described  by  Thurner  et  al. 
46|  was  sufficiently  small  (0.917  for  m  =  4, 16,  and  256)  that  we  abandoned  this  construct. 


Four  of  the  techniques  we  have  discussed  in  this  Chapter  (spectral,  wavelet,  detrended 
fluctuation  analysis,  and  Allan  factor)  yield  both  scale-independent  and  scale-dependent 
measures  and  therefore  afford  us  the  opportunity  of  directly  comparing  these  two  classes 
of  measures  in  individual  calculations:  aw  crwav(32);  as  «-»■  5V (1/32);  ^  DFA (32); 

a  a  A(10).  In  each  of  these  cases  the  fixed-scale  measure  is  found  to  greatly  outperform 
the  fractal-exponent  measure  for  all  data  sizes  examined,  as  we  previously  illustrated  for  the 
pairs  aw  crwav(32)  and  as  5V (1/32)  [pE7|[.  These  results  were  recently  confirmed  in  a 
follow-up  Israeli-Danish  study  Jd(|.  Moreover,  in  contrast  with  the  substantial  variability 
returned  in  fractal-exponent  estimates,  results  for  the  different  scale-dependent  measures  at 
m  —  32  intervals  bear  reasonable  similarity  to  each  other. 


Nunes  Amaral  et  al.  0  recently  concluded  exactly  the  opposite,  namely  that  scaling 


exponents  provide  superior  performance  to  scale-dependent  measures.  This  may  be  because 
they  relied  exclusively  on  the  distribution-dependent  measures  7]  =  h2  and  d2  (see  subsec¬ 
tions  [2.7. 1|  and  [3.5. 1|)  rather  than  on  distribution-independent  ROC  analysis.  These  same 
authors  [[56]  also  purport  to  glean  information  from  higher  moments  of  the  wavelet  coef¬ 


ficients,  but  the  reliability  of  such  information  is  questionable  because  estimator  variance 
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increases  with  moment  order  0. 


3.6.2  Computation  Times  of  the  Various  Measures 


The  computation  times  for  the  16  measures  considered  in  this  Chapter  are  provided  in 
Table  2.  All  measures  were  run  ten  times  and  averaged  except  for  the  two  DFA  measures 
which,  because  of  their  long  execution  time,  were  run  only  once.  These  long  execution  times 
are  associated  with  the  suggested  method  for  computing  DFA  |hj ,  which  is  an  N2  process. 
DFA  computation  times  therefore  increase  as  the  square  of  the  number  of  intervals  whereas 
all  14  other  methods,  in  contrast,  are  either  of  order  N  or  Nlog(N). 


Based  on  computation  time,  we  can  rank  order  the  six  scale-dependent  measures  that 
fall  into  the  first  class,  from  fastest  to  slowest:  crwav(32),  S'r(l/32),  A(10),  LF,  VLF,  and 
DFA(32).  Because  of  its  computational  simplicity  the  wavelet-transform  standard  deviation 
c wav (32)  computes  more  rapidly  than  any  of  the  other  measures.  It  is  3  times  faster  than 
that  of  its  nearest  competitor  SV(l/32),  16.5  times  faster  than  LF,  and  32500  times  faster 
than  DFA (32). 


3.6.3  Comparing  the  Most  Effective  Measures 


In  subsections  [L3|  |3.5|.  we  established  that  six  measures  succeeded  in  completely  separating 
the  normal  subjects  from  the  CHF  patients  in  our  database  for  data  lengths  L  >  215  =  32768 
RR  intervals:  VLF,  LF,  A(10),  crwav(32),  5V(l/32),  and  DFA(32).  We  now  demonstrate  from 
a  fundamental  point  of  view  that  these  measures  can  all  be  viewed  as  transformations  of  the 
interval-based  power  spectral  density  and  that  all  are  therefore  closely  related. 


Consider  first  the  relationships  of  VLF  and  LF  to  SV(l/32).  At  a  frequency  /  =  1/32  « 
0.031  cycles/interval,  Sr(f)  separates  CHF  patients  from  normal  subjects;  however,  a  range 
of  values  of  /  provide  varying  degrees  of  separability.  For  the  data  sets  we  analyzed,  sep¬ 
aration  extends  from  /  =  0.02  to  /  =  0.07  cycles/interval  [[F|.  Recall  that  VLF  and  LF 
are  simply  integrals  of  ST(f )  over  particular  ranges  of  /:  from  /  =  0.003  to  /  =  0.04  cy¬ 
cles/interval  for  VLF,  and  from  /  =  0.04  to  /  =  0.15  cycles/interval  for  LF.  Since  these 
ranges  overlap  with  those  for  which  the  power  spectral  density  itself  provides  separability,  it 
is  not  surprising  that  these  spectral  integrals  also  exhibit  this  property.  This  illustrates  the 
close  relationship  among  VLF,  LF,  and  SV(l/32). 


We  next  turn  to  the  relation  between  awav(32)  and  5V (1/32).  In  Sec.  |2.5.3|  the  relationship 
between  these  two  measures  was  established  as 
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-OO 


(23) 
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revealing  that  cr^av(m)  is  simply  a  bandpass-filtered  version  of  the  interval-based  power 
spectral  density  for  stationary  signals.  For  the  Haar  wavelet,  the  bandpass  filter  H(mf ) 
has  a  center  frequency  near  1/m,  and  a  bandwidth  about  that  center  frequency  also  near 
1/m  [[49].  Accordingly,  cr^av (32)  is  equivalent  to  a  weighted  integral  of  ST(f)  centered  about 
/  =  1/32  as  0.031,  with  an  approximate  range  of  integration  of  0.031±0.016  (0.016  to  0.047). 
Thus  the  separability  of  (x^av( 32)  is  expected  to  resemble  that  of  SV(l/32),  and  therefore  that 
of  VLF  and  LF  as  well. 


Now  consider  the  count-based  Allan  factor  evaluated  at  a  counting  time  of  ten  seconds, 
A(10).  Since  the  Allan  factor  is  directly  related  to  the  variance  of  the  Haar  wavelet  trans¬ 
form  of  the  counts  |[^J,  A(10)  is  proportional  to  the  wavelet  transform  of  the  counting 
process  evaluated  over  a  wavelet  duration  of  T  =  20  sec.  The  count-based  and  interval- 
based  measures  are  approximately  related  via  Eq.  (7)  so  that  A (10)  should  be  compared 
with  crwav(20/E[r])  ~  crwav(25)  whose  value,  in  turn,  is  determined  by  5V (1/25  =  0.04)  and 
spectral  components  at  nearby  frequencies.  Since  this  range  of  ST(f )  is  known  to  exhibit 
separability,  so  too  should  A(10). 


The  relationship  between  detrended  fluctuation  analysis  and  the  other  measures  proves 
more  complex  than  those  delineated  above.  Nevertheless,  DFA(32)  can  also  be  broadly  in¬ 
terpreted  in  terms  of  an  underlying  interval-based  power  spectral  density.  To  forge  this  link, 
we  revisit  the  calculation  of  F{m )  (see  Sec.  |2.5.4|).  The  mean  is  first  subtracted  from  the 
sequence  of  RR  intervals  {t,}  resulting  in  a  new  sequence  that  has  zero  mean.  This  new  ran¬ 
dom  process  has  a  power  spectral  density  everywhere  equal  to  ST(f),  except  at  /  =  0  where 
it  assumes  a  value  of  zero.  The  summation  of  this  zero-mean  sequence  generates  the  series 
y(k),  as  shown  in  Eq.  (0)-  Since  summation  is  represented  by  a  transfer  function  l/j2nf 
for  small  frequencies  /,  with  j  =  — 1,  the  power  spectral  density  of  y(k)  is  approximately 

given  by 


Sy(f) 


0 

Sr(f) 


f  =  0 
/  ±  0. 


(24) 


Next,  the  sequence  y(k)  is  divided  into  segments  of  length  m.  For  the  first  segment,  k  —  1 
to  m,  the  best  fitting  linear  trend  ym{k )  =  (3  +  jk  is  obtained;  the  residuals  y{k)  —  ym(k ) 
therefore  have  the  minimum  variance  over  all  choices  of  (3  and  7.  Over  the  duration  of  the 
entire  sequence  y(k),  the  local  trend  function  ym(k )  will  consist  of  a  series  of  such  linear 
segments,  and  is  therefore  piecewise  linear.  Since  ym(k )  changes  behavior  over  a  time  scale 
of  m  intervals,  its  power  will  be  concentrated  at  frequencies  below  1/m.  Since  the  residuals 
y{k )  —  ym{k )  have  minimum  variance,  the  spectrum  of  ym{k)  will  resemble  that  of  y(k)  as 
much  as  possible,  and  in  particular  at  frequencies  below  1/m  where  most  of  its  power  is 
concentrated.  Thus 


Sym(f) 


Sy(f )  I/I  <  V™ 
0  [/I  >  1/m. 


(25) 
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Equation  (^)  shows  the  mean-square  fluctuation  F2{m )  to  be  an  estimate  of  the  variance 
of  the  residuals  y(k)  —  ym(k).  By  Parseval’s  theorem,  this  variance  equals  the  integral  of  the 
spectrum  of  the  residuals  over  frequency.  Therefore 

POO  POO 

F2(m)  «  2  /  [S„(/)  -  Sam(f)]df  =  (2JT2)-1  /  Sr(f)r2df  (26) 

J 0+  Jl/m 

so  that  the  mean-square  fluctuation  F2{m)  can  be  represented  as  a  weighted  integral  of  the 
RR-interval  spectrum  over  the  appropriate  frequency  range. 

We  conclude  that  all  six  measures  that  best  distinguish  CHF  patients  from  normal  sub¬ 
jects  are  related  in  terms  of  transformations  of  the  interval-based  power  spectral  density. 


3.6.4  Selecting  the  Best  Measures 

In  the  previous  section  we  showed  that  the  six  most  effective  measures  are  interrelated  from 
a  fundamental  point  of  view.  This  conclusion  is  confirmed  by  the  overall  similarity  of  the 
ROC-area  curves  for  these  measures,  as  can  be  seen  in  Fig.  8. 

We  now  proceed  to  compare  these  ROC  areas  more  carefully  in  Fig.  9.  This  figure 
presents  1— ROC  area  plotted  against  the  data  length  analyzed  L  on  a  doubly  logarithmic 
plot.  This  is  in  contrast  to  the  usual  presentation  (such  as  Fig.  8)  in  which  ROC  area  is 
plotted  against  L  on  a  semilogarithmic  plot.  The  unconventional  form  of  the  plot  in  Fig.  9 
is  designed  to  allow  small  deviations  from  unity  area  to  be  highlighted.  Clearly,  small  values 
of  the  quantity  1— ROC  area  are  desirable.  Of  course,  as  the  data  length  L  diminishes,  ROC 
area  decreases  (and  therefore  1— ROC  area  increases)  for  all  six  measures. 

The  trends  exhibited  by  all  six  measures  displayed  in  Fig.  9  (VLF,  24(10),  DFA(32), 
crwav(32),  5V(l/32),  and  LF)  are  indeed  broadly  similar,  as  we  observed  in  Fig.  8.  For 
data  sets  that  are  sufficiently  long,  the  six  measures  provide  identical  diagnostic  capabilities 
(unity  ROC  area).  At  shorter  data  lengths,  however,  VLF,  A(10),  and  DFA(32)  exhibit  more 
degradation  than  do  the  other  three  measures.  In  accordance  with  the  analysis  provided  in 
Sec.  |3.6.3j.  VLF  deviates  the  most  at  small  values  of  L  because  it  incorporates  information 
from  values  of  ST(f )  well  below  the  frequency  range  of  separability.  Similarly  DFA(32)  favors 
low-frequency  contributions  by  virtue  of  the  l//2  weighting  function  in  Eq.  (^6j).  Finally, 
the  performance  of  A(10)  is  suppressed  because  of  the  confounding  effect  of  values  of  ST(f) 
for  /  >  0.07,  as  well  as  from  possible  errors  deriving  from  our  assumptions  regarding  the 
equivalence  of  the  count-based  and  interval-based  wavelet-transform  variances. 

The  remaining  three  measures  therefore  emerge  as  superior:  LF,  <jwav(32),  and  SV(l/32). 
Judged  on  the  basis  of  both  performance  (Fig.  9)  and  computation  time  (Table  2),  we 
conclude  that  two  of  these  are  best:  <rwav  (32)  and  S  (1/32). 
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It  is  desirable  to  confirm  these  conclusions  for  other  records,  comprising  out-of-sample 
data  sets  from  CHF  patients  and  normal  subjects.  It  will  also  be  of  interest  to  examine  the 
correlation  of  ROC  area  with  severity  of  cardiac  dysfunction  as  judged,  for  example,  by  the 
NYHA  clinical  scale.  It  is  also  important  to  conduct  comparisons  with  other  CHF  studies 
that  make  use  of  HRV  measures  [[3J . 

A  sufficiently  large  database  of  normal  and  CHF  records  would  make  it  possible  to  ex¬ 
amine  the  detailed  distinctions  between  the  two  classes  of  patients  over  a  range  of  scales 
surrounding  m  —  32.  A  superior  measure  which  draws  on  these  distinctions  in  an  optimal 
way  could  then  be  designed;  for  example  the  weighting  function  in  an  integral  over  ST(f) 
could  be  customized.  Other  measures  that  simultaneously  incorporate  properties  inherent 
in  collections  of  the  individual  measures  considered  here  could  also  be  developed. 


4  Markers  for  Other  Cardiac  Pathologies 


Returning  to  the  individual-value  plots  in  Fig.  5,  we  briefly  examine  the  behavior  of  the 
16  measures  for  several  other  cardiovascular  disorders  for  which  these  measures  may  hold 
promise.  Among  these  are  three  CHF  patients  who  also  suffer  from  atrial  fibrillation  (CHF 
c  AF,  third  column,  A),  one  heart-transplant  patient  (TRANSPLANT,  fourth  column,  $£■), 
one  ventricular-tachycardia  patient  (VT,  fifth  column,  <0),  one  sudden-cardiac-death  patient 
(SCD,  sixth  column,  □),  and  two  sleep-apnea  patients  (APNEA,  last  column,  +).  A  summary 
of  the  more  salient  statistics  of  the  original  RR  recordings  for  these  patients  is  presented  in 
Table  Al.  Prior  to  analysis,  the  RR  records  were  truncated  at  L max  =  75821  RR  intervals 
to  match  the  lengths  of  the  CHF-patient  and  normal-subject  records.  The  two  sleep-apnea 
data  records  were  drawn  from  the  MIT-BIH  Polysomnographic  Database  (Harvard-MIT 
Division  of  Health  Sciences  and  Technology)  and  comprised  16874  and  15751  RR  intervals 
respectively.  Because  of  the  limited  duration  of  the  sleep  period,  records  containing  larger 
numbers  of  RR  intervals  are  not  available  for  this  disorder. 

Too  few  patient  recordings  are  available  to  us  to  reasonably  evaluate  the  potential  of 
these  HRV  measures  for  the  diagnosis  of  these  forms  of  cardiovascular  disease,  but  the 
results  presented  here  might  suggest  which  measures  would  prove  useful  given  large  numbers 
of  records. 

We  first  note  that  for  most  of  the  16  measures,  the  values  for  the  three  CHF  patients 
with  atrial  fibrillation  (A)  tend  to  fall  within  the  range  set  by  the  other  twelve  CHF  patients 
without  atrial  fibrillation.  In  particular  the  same  six  measures  that  completely  separate 
normal  subjects  from  CHF  patients  s  AF  continue  to  do  so  when  CHF  patients  c  AF  are 
included.  The  inclusion  of  the  AF  patients  reduces  the  separability  only  for  three  of  the  16 
measures:  pNN50,  HF,  and  aw- 
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Results  for  the  heart  transplant  patient  (^)  appear  towards  the  pathological  end  of  the 
CHF  range  for  many  measures,  and  for  some  measures  the  transplant  value  extends  beyond 
the  range  of  all  of  the  CHF  patients.  Conversely,  the  sleep  apnea  values  (+)  typically 
fall  at  the  end  of,  or  beyond,  the  range  spanned  by  the  normal  subjects.  Indeed,  three 
thresholds  can  be  chosen  for  SV(l/32)  (i.e. ,  0.000008,  0.0006,  and  0.00436)  which  completely 
separate  four  classes  of  patients:  sleep  apnea  [largest  values  of  SV(l/32)],  normal  subjects 
(next  largest),  CHF  with  or  without  atrial  fibrillation  (next  largest),  and  heart  transplant 
(smallest).  While  such  separation  will  no  doubt  disappear  as  large  numbers  of  additional 
patient  records  are  included,  we  nevertheless  expect  the  largest  variability  (measured  at  32 
intervals  as  well  as  at  other  scales)  to  be  associated  with  sleep-apnea  patients  since  their 
fluctuating  breathing  patterns  induce  greater  fluctuations  in  the  heart  rate.  Conversely, 
heart-transplant  patients  (especially  in  the  absence  of  reinnervation)  are  expected  to  exhibit 
the  lowest  values  of  variability  since  the  autonomic  nervous  system,  which  contributes  to  the 
modulation  of  heart  rate  at  these  time  scales,  is  subfunctional. 

Results  for  the  ventricular  tachycardia  patient  (<})  lie  within  the  normal  range  for  most 
measures,  though  the  value  for  LF /HF  is  well  beyond  the  pathological  end  of  the  range  for 
CHF  patients.  Interestingly,  four  measures  [VLF,  HF,  SDANN,  and  A(10)]  yield  the  largest 
values  for  the  ventricular  tachycardia  patient  of  the  32  data  sets  examined.  Perhaps  A(10), 
which  exhibits  the  largest  separation  between  the  value  for  this  patient  and  those  for  the 
other  31  patients,  will  prove  useful  in  detecting  ventricular  tachycardia. 

Finally  we  consider  results  for  the  sudden-cardiac-death  patient  (□),  which  lie  between 
the  ranges  of  normal  subjects  and  CHF  patients  for  all  but  one  (LF)  of  the  six  measures 
which  separate  these  two  patient  classes.  The  pNN50  and  HF  results,  on  the  other  hand, 
reveal  values  beyond  the  normal  set,  thereby  indicating  the  presence  of  greater  variability 
on  a  short  time  scale  than  for  other  patients.  The  value  of  SDANN  for  the  SCD  patient  is 
greater  than  that  for  any  of  the  other  31  patients  so  it  is  possible  that  SDANN  could  prove 
useful  in  predicting  sudden  cardiac  death. 


5  Does  Deterministic  Chaos  Play  a  Role  in  Heart  Rate 
Variability? 


The  question  of  whether  normal  HRV  arises  from  a  low-dimensional  attractor  associated 
with  a  deterministic  dynamical  system  or  has  stochastic  underpinnings  has  continued  to 
entice  researchers  |H]].  Previous  studies  devoted  to  this  question  have  come  to  conflicting 
conclusions.  The  notion  that  the  normal  human-heartbeat  sequence  is  chaotic  appears  to 


have  been  first  set  forth  by  Goldberger  and  West 
as  well  as  in  opposition  to  it  [| 


Papers  in  support  of  this  hypothesis 
7T| ,  have  appeared  recently.  An  important  factor  that 


pertains  to  this  issue  is  the  recognition  that  correlation,  rather  than  nonlinear  dynamics,  can 
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lie  behind  these  seemingly  chaotic  recordings  || . 

We  address  this  question  in  this  section  and  conclude,  for  both  normal  subjects  and 
patients  with  several  forms  of  cardiac  dysfunction,  that  the  sequence  of  heartbeat  intervals 
does  not  exhibit  chaos. 


5.1  Methods 

5.1.1  Phase-Space  Reconstruction 


One  way  of  approaching  the  question  of  chaos  in  heartbeat  recordings  is  in  terms  of  a 
phase-space  reconstruction  and  the  use  of  an  algorithm  to  estimate  one  or  more  of  its  fractal 
dimensions  Dq  ||  [59],  [72|] .  The  box-counting  algorithm  [[73|]  provides  a  method  for  estimating 
the  capacity  (or  box-counting)  dimension  D0  of  the  attractor  ||  [7^,  [74|] . 


The  phase-space  reconstruction  approach  operates  on  the  basis  that  the  topological  prop¬ 
erties  of  the  attractor  for  a  dynamical  system  can  be  determined  from  the  time  series  of  a 

single  observable  [[75],  [76],  [77).  A  p-dimensional  vector  T  =  {r*,  Ti+i,  ...Tj+(p_i)z}  is  formed 
from  the  sequence  of  RR  intervals  {t*}.  The  parameter  p  is  the  embedding  dimension,  and 
l  is  the  lag,  usually  taken  to  be  the  location  of  the  first  zero  crossing  of  the  autocorrelation 


function  of  the  time  series  under  study.  As  the  time  index  i  progresses,  the  vector  T  traces 
out  a  trajectory  in  the  p-dimensional  embedding  space. 


The  box-counting  algorithm  estimates  the  capacity  dimension  of  this  trajectory.  Specif¬ 
ically,  the  negative  slope  of  the  logarithm  of  the  number  of  full  boxes  versus  the  logarithm 
of  the  width  of  the  boxes,  over  a  range  of  4/[max(r)  —  min(r)]  to  32/[max(r)  —  min(r)] 
boxes,  in  powers  of  two,  provides  an  estimate  of  Do-  For  convenience  in  the  generation  of 
phase-randomized  surrogates,  data  sets  were  limited  to  powers  of  two:  21(>  =  65536  intervals 
for  all  but  the  two  sleep  apnea  data  sets,  for  which  21 1  =  16384  and  213  =  8192  intervals  were 
used  respectively  (the  sleep  apnea  data  sets  are  much  shorter  than  the  others  as  indicated  in 
Table  Al).  For  uncorrelated  noise,  the  capacity  dimension  Do  continues  to  increase  as  the 
embedding  dimension  p  increases.  For  an  attractor  in  a  deterministic  system,  in  contrast,  Do 
saturates  as  p  becomes  larger  than  2 D0  + 1  |74J .  Such  saturation,  however,  is  not  a  definitive 
signature  of  deterministic  dynamics;  temporal  correlation  in  a  stochastic  process  can  also  be 
responsible  for  what  appears  to  be  underlying  deterministic  dynamics  [|[  [78].  [79],  |(J.  The 
distribution  of  the  data  can  also  play  a  role. 


Surrogate  data  analysis,  discussed  in  Sec.  |5.1.3|,  is  useful  in  establishing  the  underlying 
cause,  by  proving  or  disproving  various  null  hypotheses  that  the  measured  behavior  arises 
from  other,  non-chaotic  causes. 


35 


5.1.2  Removing  Correlations  in  the  Data 


Adjacent  intervals  in  heartbeat  sequences  prove  to  exhibit  large  degrees  of  correlation.  The 
serial  correlation  coefficient  for  RR  intervals 

=  E[Ti+iTj]  -  E2 \Tj\ 

P[  ’  ~  E[rf]  -  E2[tj]  ’ 

exceeds  0.9  for  two  thirds  of  the  data  sets  examined,  with  an  average  value  of  0.83;  the 
theoretical  maximum  is  unity,  with  zero  representing  a  lack  of  serial  correlation.  Since 
such  large  correlation  is  known  to  interfere  with  the  detection  of  deterministic  dynamics, 
we  instead  employ  the  first  difference  of  the  RR  intervals:  n*  =  r*  —  rl+\ .  Such  a  simple 
transformation  does  not  change  the  topological  properties  of  dynamical  system  behavior, 
and  a  related  but  more  involved  procedure  is  known  to  reduce  error  in  estimating  other 
fractal  dimensions  [|3T|. 

The  serial  correlation  coefficient  for  the  differenced  sequence  {u*}  turns  out  to  have  a 
mean  value  of  -0.084,  and  none  exceeds  0.6  in  magnitude.  This  sequence  will  therefore  be 

used  to  generate  the  embedding  vector  T  =  {n*,  v*+i,  ...nj+p_i},  with  the  lag  l  set  to  unity 
since  p{v)  ~  0. 


5.1.3  Surrogate  Data  Analysis 

Information  about  the  nature  of  an  interval  or  segment  series  may  be  obtained  by  applying 
various  statistical  measures  to  surrogate  data  sets.  These  are  processes  constructed  from  the 
original  data  in  ways  designed  to  preserve  certain  characteristics  of  the  original  data  while 
eliminating  (or  modifying)  others.  Surrogate  data  analysis  provides  a  way  of  determining 
whether  a  given  result  arises  from  a  particular  property  of  the  data  set. 

We  make  use  of  two  kinds  of  surrogate  data  sets:  shuffled  intervals  and  randomized 
phases.  In  particular,  we  compare  statistical  measures  calculated  from  both  the  original 
data  and  from  its  surrogates  to  distinguish  those  properties  of  the  data  set  that  arise  from 
correlations  (such  as  from  long-term  rate  fluctuations)  from  those  properties  inherent  in  the 
form  of  the  original  data. 


Shuffled  intervals.  The  first  class  of  surrogate  data  is  formed  by  shuffling  (randomly 
reordering)  the  sequence  of  RR  intervals  of  the  original  data  set.  Such  random  reordering 
destroys  dependencies  among  the  intervals,  and  therefore  the  correlation  properties  of  the 
data,  while  exactly  preserving  the  interval  histogram.  As  a  refinement,  it  is  possible  to 
divide  the  data  into  blocks  and  shuffle  intervals  within  each  block,  while  keeping  each  block 
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separate.  In  this  case  dependencies  remain  over  durations  much  larger  than  the  block  size, 
while  being  destroyed  for  durations  smaller  than  it. 


Randomized  phases.  The  other  class  of  surrogate  data  we  consider  is  obtained  by  Fourier 
transforming  the  original  interval  data,  and  then  randomizing  the  phases  while  leaving  the 
spectral  magnitudes  intact.  The  modified  function  is  then  inverse-transformed  to  return  to 
a  time-domain  representation  of  the  intervals.  This  technique  exactly  preserves  the  second- 
order  correlation  properties  of  the  interval  sequence  while  removing  other  temporal  structure, 
for  example  that  needed  for  phase-space  reconstruction.  The  interval  histogram  typically 
becomes  Gaussian  for  this  surrogate  as  a  result  of  the  central  limit  theorem. 


5.2  Absence  of  Chaos 

In  Fig.  10  we  present  plots  of  the  capacity  dimension  D0  plotted  against  the  embedding 
dimension  p  for  the  12  normal  subjects.  Similar  plots  are  presented  in  Fig.  11  for  CHF 
patients  without  atrial  fibrillation,  and  in  Fig.  12  for  patients  with  other  cardiac  pathologies. 
For  most  of  the  panels  in  Fig.  10,  the  D0  estimate  for  the  original  data  (solid  curve)  is  indeed 
seen  to  rise  with  increasing  embedding  dimension  p,  and  shows  evidence  of  approaching  an 
asymptotic  value  of  about  2.  The  phase-randomized  data  (dotted  curve)  yields  somewhat 
higher  results,  consistent  with  the  presence  of  a  deterministic  attractor. 

However,  estimates  of  D0  for  the  shuffled  data  (dashed  curve)  nearly  coincide  with  those 
of  the  original  data.  This  surrogate  precisely  maintains  the  distribution  of  the  relative  sizes 
of  the  differenced  RR  intervals,  but  destroys  any  correlation  or  other  dependencies  in  the 
data.  Specifically,  no  deterministic  structure  remains  in  the  shuffled  surrogate.  Therefore, 
the  apparent  asymptote  in  Dq  must  derive  from  the  particular  form  of  the  RR-interval  dis¬ 
tribution,  and  not  from  a  nonlinear  dynamical  attractor.  To  quantify  the  closeness  between 
results  for  original  and  shuffled  data,  the  shuffling  process  was  performed  ten  times  with  ten 
different  random  seeds.  The  average  of  these  ten  independent  shufflings  forms  the  dashed 
curve,  with  error  bars  representing  the  largest  and  smallest  estimates  obtained  for  Do  at 
each  corresponding  embedding  dimension.  Results  for  the  original  data  indeed  lie  within 
the  error  bars.  Therefore,  the  original  data  does  not  differ  from  the  shuffled  data  at  the 
[(10  —  1)/(10  +  1)]  100%  =  82%  significance  level. 

Results  for  a  good  number  of  the  other  31  data  sets,  normal  and  pathological  alike,  are 
nearly  the  same;  the  original  data  lie  within  the  error  bars  of  the  shuffled  surrogates  and 
consequently  do  not  exhibit  significant  evidence  of  deterministic  dynamics.  However  for 
some  records  the  original  data  exceed  the  limits  set  by  the  ten  shufflings.  But  attributing 
this  behavior  to  low-dimensional  chaos  proves  difficult  for  almost  all  the  data  records,  for 
three  reasons.  First,  the  surrogates  often  yield  smaller  results  for  Dq  than  the  original  data, 
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while  shuffling  data  from  a  true  chaotic  system  will  theoretically  yield  larger  Dq  estimates. 
Second,  of  those  data  sets  for  which  the  original  results  fall  below  those  of  their  shuffled 
surrogates,  few  achieve  an  asymptotic  value  for  Do,  but  rather  continue  to  climb  steadily 
with  increasing  embedding  dimension.  Finally,  for  the  remaining  data  sets,  the  original  and 
shuffled  data  yield  curves  that  do  not  differ  nearly  as  much  as  do  those  from  simulated  systems 
known  to  be  chaotic,  such  as  the  Henon  attractor  (not  shown)  [^2j.  Taken  together,  these 
observations  suggest  that  finite  data  length,  residual  correlations,  and/or  other  unknown 
effects  in  these  RR-interval  recordings  give  rise  to  behavior  in  the  phase-space  reconstruction 
that  can  masquerade  as  chaos. 

Results  for  the  phase-randomized  surrogates  (dotted  curves  in  Figs.  10-12),  while  exhibit¬ 
ing  a  putative  capacity  dimension  in  excess  of  the  original  data  and  of  the  shuffled  surro¬ 
gates  for  every  patient,  exhibit  a  remarkable  similarity  to  each  other.  This  occurs  because 
the  differenced-interval  sequence  {vi}  displays  little  correlation  and  the  phase-randomization 
process  preserves  this  lack  of  correlation  while  imposing  a  Gaussian  differenced-interval  dis¬ 
tribution.  The  result  is  nearly  white  Gaussian  noise  regardless  of  the  original  data  set,  and 
indeed  results  for  independently  generated  white  Gaussian  noise  closely  follow  these  curves 
(not  shown). 


Time  series  measured  from  the  heart  under  distinctly  non-normal  operating  conditions 
such  as  ventricular  fibrillation  [S3  ,  or  those  obtained  under  non-physiological  conditions  such 
as  from  excised  hearts  |p4| ,  may  exhibit  chaotic  behavior  under  some  circumstances.  This  is 


precisely  the  situation  for  cellular  vibrations  in  the  mammalian  cochlea  |55|  [33].  [sT|| .  When 


the  exciting  sound  pressure  level  (SPL)  is  in  the  normal  physiological  regime,  the  cellular- 
vibration  velocities  are  nonlinear  but  not  chaotic.  When  the  SPL  is  increased  beyond  the 
normal  range,  however,  routes  to  chaos  emerge.  It  has  been  suggested  that  since  chaos- 
control  methods  prove  effective  for  removing  such  behavior,  the  heartbeat  sequence  must  be 
chaotic.  This  is  an  unwarranted  conclusion  since  such  methods  often  work  equally  well  for 
purely  stochastic  systems,  which  cannot  be  chaotic 


6  Mathematical  Models  for  Heart  Rate  Variability 


The  emphasis  in  this  Chapter  has  been  on  HRV  analysis.  The  diagnostic  capabilities  of 
various  measures  have  been  presented  and  compared.  However  the  results  that  emerge  from 
these  studies  also  serve  to  expose  the  mathematical-statistical  character  of  the  RR  time 
series.  Such  information  can  be  of  benefit  in  the  development  of  mathematical  models  which 
can  be  used  to  simulate  RR  sequences.  Such  models  may  prove  useful  for  generating  realistic 
heartbeat  sequences  in  pacemakers,  for  example,  and  for  developing  improved  physiological 
models  of  the  cardiovascular  system.  In  this  section  we  evaluate  the  suitability  of  several 
integrate-and-fire  constructs  for  modeling  heart  rate  variability. 
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6.1  Integrate-and-Fire  Model 


Integrate-and-fire  models  are  widely  used  in  the  neurosciences  [|9|  as  well  as  in  cardiology 
IP  0.  [n].  []2|.  |93| .  These  models  are  attractive  in  part  because  they  capture  known  physiology 
in  a  simple  way.  The  integration  of  a  rate  process  can  represent  the  cumulative  effect  of 
neurotransmitter  on  the  postsynaptic  membrane  of  a  neuron,  or  the  currents  responsible 
for  the  pacemaker  potential  in  the  sino-atrial  node  of  the  heart.  The  crossing  of  a  preset 
threshold  by  the  integrated  rate  then  gives  rise  to  an  action  potential,  or  a  heart  contraction. 

The  sequence  of  discrete  events  comprising  the  human  heartbeat  can  be  viewed  as  a  point 
process  deriving  from  a  continuous  rate  function.  The  integrate-and-fire  method  is  perhaps 
the  simplest  means  for  generating  a  point  process  from  a  rate  process  [S.  17].  In  this  model, 


illustrated  schematically  in  Fig.  13,  the  rate  function  A (t)  is  integrated  until  it  reaches  a  fixed 
threshold  9,  whereupon  a  point  event  is  generated  and  the  integrator  is  reset  to  zero.  Thus 
the  occurrence  time  of  the  ( i  +  l)st  event  is  implicitly  obtained  from  the  first  occurrence  of 


Ei+1 


A (t)  dt  =  9. 


(27) 


The  mean  heart  rate  A  (beats/sec)  is  given  by  A  =  E[A(t)]  =  1/E[r].  Because  the  conversion 
from  the  rate  process  to  the  point  process  is  so  direct,  most  statistics  of  the  point  process 
closely  mimic  those  of  the  underlying  rate  process  A (f).  In  particular,  for  frequencies  sub¬ 
stantially  lower  than  the  mean  heart  rate,  the  theoretical  power  spectral  densities  of  the 
point  process  and  the  underlying  rate  coincide 


17 


6.2  Kernel  of  the  Integrate-and-Fire  Model 

The  detailed  statistical  behavior  of  the  point  process  generated  by  the  integrate-and-fire 
construct  depends  on  the  statistical  character  of  the  rate  function  in  the  integrand  of  Eq.  (|7j) 
p.  The  fractal  nature  of  the  heartbeat  sequence  requires  that  the  rate  function  A (t)  itself 
be  fractal.  Assuming  that  the  underlying  rate  is  stochastic  noise,  there  are  three  plausible 
candidates  for  \(t)  Q:  fractal  Gaussian  noise,  fractal  lognormal  noise,  and  fractal  binomial 
noise.  We  briefly  discuss  these  three  forms  of  noise  in  turn. 


6.2.1  Fractal  Gaussian  Noise 

Fractal  Gaussian  noise  served  admirably  in  our  initial  efforts  to  develop  a  suitable  integrate- 
and-fire  model  for  the  heartbeat  sequence  p.  A  Gaussian  process  has  values  that,  at  any 
set  of  times,  form  a  jointly  Gaussian  random  vector.  This  property,  along  with  its  mean  and 
spectrum,  completely  define  the  process.  We  consider  a  stationary  rate  process  with  a  mean 
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equal  to  the  expected  heart  rate,  and  a  l//-type  rate-based  spectrum  that  takes  the  form 

E3 

=  A  2i(/tlme)  +  A[1  +  (/tlm„//o)-“].  (28) 


Here,  A  (beats/sec)  and  a  are  the  previously  defined  heart  rate  and  fractal  exponent  re¬ 
spectively,  <5(-)  is  again  the  Dirac  delta  function,  and  /o  is  the  cutoff  frequency  (cycles/sec). 
Mathematically,  cutoffs  at  low  or  high  frequencies  (or  sometimes  both)  are  required  to  ensure 
that  r  has  a  finite  variance.  In  practice,  the  finite  duration  of  the  simulation  guarantees  the 
former  cutoff,  while  the  finite  time  resolution  of  the  simulated  samples  of  r  guarantees  the 
latter.  It  is  to  be  noted  that  the  designation  fractal  Gaussian  noise  properly  applies  only  for 
a  <  1;  the  range  1  <  a  <  3  is  generated  by  fractal  Brownian  motion  [19).  In  the  interest  of 
simplicity,  however,  we  employ  the  term  fractal  Gaussian  noise  for  all  values  of  a. 


There  are  a  number  of  recipes  in  the  literature  for  generating  fractal  Gaussian  noise 
,  |97]|.  All  typically  result  in  a  sampled  version  of  fractal  Gaussian  noise,  with 
equal  spacing  between  samples.  A  continuous  version  of  the  signal  is  obtained  by  using  a 
zeroth-order  interpolation  between  the  samples. 


The  use  of  a  fractal  Gaussian  noise  kernel  in  the  integrate-and-hre  construct  results  in 
the  so-called  fractal-Gaussian-noise  driven  integrate-and-hre  model  [H  |17|].  The  mean  of  the 


rate  process  is  generally  required  to  be  much  larger  than  its  standard  deviation,  in  part  so 
that  the  times  when  the  rate  is  negative  (during  which  no  events  may  be  generated)  remain 
small. 


The  interval- based  spectrum  is  obtained  by  applying  Eq.  (0)  to  Eq.  (|8D  (see  Sec.  ^[2] 
and  Ref.  [:9J).  The  following  approximate  result,  in  terms  of  interval-based  frequency  /,  is 
obtained: 


st(/)«a-W)  +  i  +  (a///„)-“], 


(29) 


We  proceed  to  discuss  two  other  physiologically  based  noise  inputs 
fractal  Gaussian  noise  in  appropriate  limits. 


17  ;  both  reduce  to 


6.2.2  Fractal  Lognormal  Noise 

A  related  process  results  from  passing  fractal  Gaussian  noise  through  a  memoryless  ex¬ 
ponential  transform.  This  process  plays  a  role  in  neurotransmitter  exocytosis  |)8|].  Since 
the  exponential  of  a  Gaussian  is  lognormal,  we  refer  to  this  process  as  fractal  lognormal 
noise  0,  El-  If  X(t)  denotes  a  fractal  Gaussian  noise  process  with  mean  E[X],  variance 
Var[AT],  and  autocovariance  function  Kx(t),  then  \(t)  =  exp[X(f)]  is  a  fractal  lognormal 
noise  process  with  moments  E[An]  =  exp(?rE[A"]  +  n2Var[A"]/2)  and  autocovariance  function 
Kx(t)  =  E2 [A] (exp[A\- (r )]  -  1)  [0. 
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By  virtue  of  the  exponential  transform,  the  autocorrelation  functions  of  the  Gaussian 
process  X(t)  and  the  lognormal  process  A (t)  differ;  thus  their  spectra  differ.  In  particular, 
the  power  spectral  density  of  the  resulting  point  process  does  not  follow  the  exact  form  of 
Eq.  (^),  although  for  small  values  of  Var[X]  the  experimental  periodogram  closely  resembles 
this  ideal  form  |5S|.  In  the  limit  of  very  small  values  of  Var[X],  the  exponential  operation 
approaches  a  linear  transform,  and  the  rate  process  reduces  to  fractal  Gaussian  noise. 


6.2.3  Fractal  Binomial  Noise 


A  third  possible  kernel  for  the  integrate-and-fire  model  is  fractal  binomial  noise.  This  process 
results  from  the  addition  of  a  number  of  independent  identical  alternating  fractal  renewal  pro¬ 
cesses.  The  sum  is  a  binomial  process  with  the  same  fractal  exponent  as  each  of  the  individual 
alternating  fractal  processes  [EM  |100|| .  This  binomial  process  can  serve  as  a  rate  function  for 


an  integrate-and-fire  process;  the  result  is  the  fractal-binomial-noise  driven  integrate-and-fire 
model  1 17].  This  construct  was  initially  designed  to  model  the  superposition  of  alternating 


currents  from  a  number  of  ion  channels  with  fractal  behavior  |]99|,  I100L II Oil,  110211 .  As  the  num¬ 


ber  of  constituent  processes  increases,  fractal  binomial  noise  converges  to  fractal  Gaussian 
noise  with  the  same  fractal  exponent  therefore  leading  to  the  fractal-Gaussian-noise  driven 
integrate-and-fire  point  process 
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6.3  Jittered  Integrate-and-Fire  Model 


Conventional  integrate-and-fire  constructs  have  only  a  single  source  of  randomness,  the  rate 
process  A (f).  It  is  sometimes  useful  to  incorporate  a  second  source  of  randomness  into 
such  models.  One  way  to  carry  this  out  is  to  impose  random  jitter  on  the  interevent  times 
generated  in  the  integrate-and-fire  process 
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The  procedure  used  for  imposing  this  jitter  is  as  follows.  After  generating  the  time  of  the 


( i  +  l)st  event  ti+ i  in  accordance  with  Eq.  (p7|),  the  ( i  +  l)st  interevent  time  Tj+i  =  ti+ 1  —  G 
is  multiplied  by  a  dimensionless  Gaussian-distributed  random  variable  with  unit  mean  and 
variance  cr2.  This  results  in  the  replacement  of  ti+ 1  by 


ti+ i  +  a(ti+ 1  -  ti)Mi  (30) 

before  the  simulation  of  subsequent  events  {ti+ 2,  fi+3,. ..).  The  quantity  (A/)}  represents 
a  sequence  of  zero-mean,  unity-variance  independent  Gaussian  random  variables  and  the 
standard  deviation  cr  is  a  free  parameter  that  controls  the  strength  of  the  proportional  jitter. 
In  accordance  with  the  construction  of  the  model,  events  at  subsequent  times  experience 
jitter  imposed  by  all  previous  events  (see  Fig.  14). 
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The  overall  result  is  the  jittered  integrate-and-fire  model.  The  jitter  serves  to  introduce  a 
source  of  scale-independent  noise  into  the  point  process.  Depending  on  the  character  of  the 
input  rate  function,  this  model  can  give  rise  to  the  fractal-Gaussian-noise,  fractal-lognormal- 
noise,  or  fractal-binomial-noise  driven  jittered  integrate-and-fire  processes,  among  others. 
Two  limiting  behaviors  readily  emerge  from  this  model.  When  a  —*  0,  the  jitter  disappears 
and  the  jittered  integrate-and-fire  model  reduces  to  the  ordinary  integrate-and-fire  model. 
At  the  other  extreme,  when  cr  — >  oo,  the  jitter  dominates  and  the  result  is  a  homogeneous 
Poisson  process.  The  fractal  behavior  present  in  the  rate  function  A (t)  then  disappears  and 
\(t)  behaves  as  if  it  were  constant.  Between  these  two  limits,  as  cr  increases  the  fractal  onset 
time  l//o  of  the  resulting  point  process  increases  and  the  fractal  characteristics  of  the  point 
process  are  progressively  lost,  first  at  higher  frequencies  (shorter  times)  and  subsequently  at 
lower  frequencies  (longer  times)  fTT|. 


6.3.1  Simulating  the  Jittered  Integrate-and-Fire  Point  Process 

We  use  a  fractal  Gaussian  noise  kernel  in  the  jittered  integrate-and-fire  construct  to  simulate 
the  human  heartbeat.  Fractal  Gaussian  noise  with  the  appropriate  mean,  standard  deviation, 
and  fractal  exponent  as  is  generated  using  Fourier-transform  methods.  This  noise  serves  as 
the  kernel  in  an  integrate-and-fire  element  [Eq.  (^7])],  the  output  of  which  is  jittered  in 
accordance  with  Eq.  (|50|).  For  each  of  the  12  normal  subjects  and  15  CHF  (s  and  c  AF) 
patients,  the  model  parameters  were  adjusted  iteratively  to  obtain  a  simulated  data  set  that 
matched  the  patient  recording  as  closely  as  possible.  In  particular,  each  of  the  27  simulations 
contained  exactly  Lmax  =  75821  RR  intervals,  displayed  mean  heart  rates  that  fell  within  1% 
of  the  corresponding  data  recordings,  and  exhibited  wavelet-transform  standard  deviations 
that  were  nearly  identical  with  those  of  the  data.  The  simulation  results  for  CHF  patients 
with  and  without  atrial  fibrillation  were  very  similar. 


6.3.2  Statistics  of  the  Simulated  Point  Process  for  Normal  Subjects  and  CHF 
Patients 

Having  developed  the  jittered  integrate-and-fire  model  and  simulated  sequences  of  RR  inter¬ 
vals,  we  proceed  to  examine  some  of  the  statistical  features  of  the  resulting  point  processes. 

Figure  15  displays  representative  results  in  which  the  data  (solid  curves)  and  simulations 
(dotted  curves)  are  compared  for  a  single  normal  subject  (left  column)  and  a  single  CHF  (s 
AF)  patient  (right  column).  Figure  15(a)  illustrates  that  the  RR-interval  sequences  obtained 
from  the  model  and  the  data  are  qualitatively  similar,  for  both  the  normal  and  heart-failure 
cases.  Furthermore,  the  agreement  between  the  simulated  and  actual  wavelet-transform 
standard-deviation  curves  is  excellent  in  both  cases,  as  is  evident  in  Fig.  15(b).  Even  the 
gentle  increase  of  <7wav  for  the  heart-failure  patient  at  small  scale  values  is  reproduced,  by 
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virtue  of  the  jitter  introduced  into  the  model.  It  is  apparent  in  Fig.  15(c),  however,  that 
the  simulated  RR-interval  histogram  for  the  heart-failure  patient  does  not  fit  the  actual 
histogram  well;  it  is  too  broad  which  indicates  that  crint  of  the  simulation  is  too  large.  Finally, 
Fig.  15(d)  shows  that  the  simulated  spectra  of  the  RR  intervals  match  the  experimental 
spectra  quite  nicely,  including  the  whitening  of  the  heart-failure  spectrum  at  high  frequencies. 
This  latter  effect  is  the  counterpart  of  the  flattening  of  crwav  at  small  scales,  and  results  from 
the  white  noise  introduced  by  the  jitter  at  high  frequencies. 


6.3.3  Simulated  Individual  Value  Plots  and  ROC- Area  Curves 


To  examine  how  well  the  jittered  integrate-and-fire  model  mimics  the  collection  of  data  from 
a  global  perspective,  we  constructed  individual-value  plots  and  ROC-area  curves  using  the 
27  simulated  data  sets.  The  results  are  presented  in  Figs.  16  and  17  respectively.  Direct 
comparison  should  be  made  with  Figs.  5  and  8,  respectively,  which  provide  the  same  plots 
for  the  actual  heart-failure  patient  and  normal-subject  data. 


Comparing  Fig.  16  with  Fig.  5,  and  Fig.  17  with  Fig.  8,  we  see  that  in  most  cases 
the  simulated  and  actual  results  are  remarkably  similar.  Based  on  their  overall  ability  to 
distinguish  between  CHF  and  normal  simulations  for  the  full  set  of  75821  RR  intervals, 
the  16  measures  portrayed  in  Fig.  16  roughly  fall  into  three  classes.  A  similar  division  was 
observed  for  the  actual  data,  as  discussed  in  Sec.  EBB 


As  is  evident  in  Fig.  16,  five  measures  (indicated  by  boldface  font)  succeed  in  fully 
separating  the  two  kinds  of  simulations  and,  by  definition,  fall  into  the  first  class:  VLF,  LF, 
<jwav(32),  SV(l/32),  and  DFA(32).  These  measures  share  a  dependence  on  a  single  scale,  or 
small  range  of  scales,  near  32  heartbeat  intervals.  However,  of  these  five  measures  only  two 
(LF  and  VLF)  successfully  separate  the  simulated  results  for  the  next  smaller  number  of  RR 
intervals  ( L  =  65536)  and  none  succeeds  in  doing  so  for  yet  smaller  values  of  L. 


These  same  five  measures  also  fully  distinguish  the  normal  subjects  from  the  CHF  patients 
(see  Fig.  5);  however  a  sixth  measure,  A(10),  also  achieves  this.  For  the  actual  data  all  six 
of  these  measures  successfully  separate  the  actual  heart-failure  patients  from  the  normal 
subjects  for  data  lengths  down  to  L  =  32768  RR  intervals. 


The  simulated  ROC-area  curves  presented  in  Fig.  17  resemble  the  curves  for  the  normal 
and  CHF  data  shown  in  Fig.  8,  but  there  are  a  few  notable  distinctions.  The  simulated  and 
actual  ROC  areas  (Figs.  17  and  8  respectively)  for  the  interval  measures  SDANN  and  (Tint 
are  just  about  the  same  for  L  =  64  RR  intervals.  However  both  of  these  measures  exhibit 
decreasing  ROC  areas  as  the  number  of  simulated  RR  intervals  ( L )  increases  (Fig.  17).  In 
contrast,  the  ROC  areas  based  on  the  actual  data  increase  with  increasing  data  length  L 
in  normal  fashion  (Fig.  8).  The  decrease  in  ROC  area  observed  in  Fig.  17  means  that 
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performance  is  degraded  as  the  simulated  data  length  increases.  This  paradoxical  result 
likely  arises  from  a  deficiency  of  the  model;  it  proved  impossible  to  simultaneously  fit  the 
wavelet-transform  standard  deviation  at  all  scales,  the  mean  heart  rate,  and  the  RR-interval 
standard  deviation.  Choosing  to  closely  fit  the  two  former  quantities  to  the  data  rendered 
inaccurate  the  simulated  interval  standard  deviation.  Since  shorter  hies  have  insufficient 
length  to  exhibit  the  longer  term  fluctuations  that  dominate  fractal  activity,  the  interval 
standard  deviation  for  small  values  of  L  is  not  influenced  by  these  fluctuations.  Apparently 
the  manner  in  which  the  long-term  effects  are  expressed  differs  in  the  original  data  and  in 
the  simulation. 

The  most  surprising  distinction  between  the  actual  and  simulated  results,  perhaps,  is  the 
dramatic  putative  improvement  in  the  ability  of  the  five  scale-independent  fractal-exponent 
measures  to  separate  simulated  heart-failure  and  normal  data,  as  the  number  of  RR  intervals 
increases.  The  Allan  factor  exponent  a  a  appears  to  offer  the  greatest  “improvement”,  as 
is  seen  by  comparing  Figs.  17  and  8.  All  of  the  exponents  except  aw  attain  an  ROC  area 
exceeding  0.97  for  Lmax  =  75821  RR  intervals.  Nevertheless,  the  improved  separation  yields 
performance  that  remains  inferior  to  that  of  the  fixed-scale  measures.  The  apparent  improved 
separation  in  the  simulated  data  may  therefore  arise  from  a  nonlinear  interaction  in  the  model 
between  the  fractal  behavior  and  the  signal  component  near  a  scale  of  32  heartbeat  intervals. 
Or  it  may  be  that  some  aspect  of  the  data,  not  reliably  detected  by  fractal-exponent  measures 
applied  to  the  original  data,  emerges  more  clearly  following  the  modeling  and  simulation. 
This  could  mean  that  another  measure  of  the  fractal  exponent  of  the  actual  data,  perhaps 
incorporating  in  some  way  the  processing  provided  by  the  modeling  procedure,  might  yield 
better  results  than  the  more  straightforward  fractal-exponent  measures  that  we  have  used 
to  this  point. 


6.3.4  Limitations  of  the  Jittered  Integrate-and-Fire  Model 

Figure  18  presents  a  direct  comparison  between  four  measures  derived  from  the  original  data 
and  from  the  corresponding  simulations.  Each  panel  contains  12  normal  results  (+),  12  CHF 
s  AF  results  (x),  and  3  CHF  c  AF  results  (A).  The  diagonal  lines  correspond  to  perfect 
agreement. 

The  results  for  SV(l/32)  and  crwav(32)  are  excellent,  illustrating  the  remarkable  ability  of 
the  model  to  mimic  the  data  at  a  time  scale  of  32  interbeat  intervals.  The  correlation  coeffi¬ 
cients  for  these  two  measures  lie  very  close  to  unity,  indicating  almost  perfect  correspondence 
between  the  original  and  simulated  data. 

The  upper  left  panel  presents  values  for  crint ,  the  standard  deviation  of  the  RR  intervals. 
Agreement  between  the  simulation  and  original  data  is  only  fair,  with  simulated  values  of 
(Tint  generally  somewhat  higher  than  desired.  The  correlation  coefficient  p  =  0.71  for  all  27 
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simulations  indicates  significant  correlation  (p  <  0.0003  that  27  random  pairs  of  numbers 
with  the  same  distribution  would  achieve  a  magnitude  of  p  >  0.71),  but  not  extremely  close 
agreement  between  simulated  and  original  RR-interval  sequences. 

This  disagreement  highlights  a  problem  with  the  simple  jittered  integrate-and-hre  model. 
The  fractal-Gaussian-noise  integrate-and-hre  model  without  jitter  yields  substantially  closer 
agreement  between  the  data  and  simulation  for  this  statistic  [|] .  Although  the  introduction 
of  jitter  leads  to  improved  agreement  for  crwav,  it  brings  the  model  results  further  away  from 
the  data  as  measured  by  crint.  Apparently,  a  method  other  than  jitter  must  be  devised  to 
forge  increased  agreement  with  awav,  while  not  degrading  the  agreement  with  crint -  One 
possibility  is  the  reordering  of  the  individual  RR  intervals  over  short  time  scales. 

The  simulated  exponent  aw  also  fails  to  show  close  agreement  with  the  data,  suggesting 
that  there  are  other  features  of  the  model  that  are  less  than  ideal.  However,  since  this 
measure  is  of  little  use  in  separating  CHF  patients  from  normal  subjects,  this  disagreement 
takes  on  reduced  importance. 


6.4  Toward  an  Improved  Model  of  Heart  Rate  Variability 


Agreement  between  the  simulations  and  data  would  most  likely  be  improved  by  adding 
other  inputs  to  the  model  aside  from  fractal  noise,  as  illustrated  schematically  in  Fig.  19. 
Input  signals  related  to  respiration  and  blood  pressure  are  two  likely  candidates  since  it  is 
well  known  that  variations  in  these  physiological  functions  influence  heart  rate  variability. 
Experience  also  teaches  us  that  there  is  a  measure  of  white  Gaussian  noise  present,  possibly 
associated  with  the  measurement  process. 


As  indicated  in  Sec.  |6.3.4|,  the  introduction  of  jitter  in  the  integrate-and-hre  construct 
improves  the  agreement  of  the  model  with  certain  features  of  the  the  data,  but  degrades  it  for 
others.  An  adaptive  reset  will  provide  a  more  flexible  alternative  to  the  jitter.  Moreover,  the 
threshold  6  could  be  converted  into  a  stochastic  process  as  a  way  of  incorporating  variability 
in  other  elements  of  the  system  [0- 


All  things  considered,  the  jittered  integrate-and-hre  construct  does  a  rather  good  job  of 
mimicking  the  actual  data,  though  the  introduction  of  a  more  physiologically  based  model 
would  be  a  welcome  addition. 
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7  Conclusion 


For  the  purposes  of  heart-rate-variability  analysis,  the  occurrence  times  of  the  sequence 
of  R  phases  in  the  electrocardiogram  can  be  represented  as  a  fractal-rate  stochastic  point 
process.  Using  a  collection  of  standard  and  novel  heart-rate- variability  measures,  we  have 
examined  the  sequence  of  times  between  the  R  phases,  viz.  the  RR  time  series,  of  ECGs 
recorded  from  normal  subjects  and  congestive  heart-failure  patients,  as  well  as  from  several 
patients  with  other  cardiovascular  disorders.  Congestive-heart-failure  patients  who  also  suf¬ 
fered  from  atrial  fibrillation  were  treated  as  a  separate  class.  We  examined  sixteen  heart-rate- 
variability  measures,  comprising  frequency-domain,  time-domain,  scaling-exponent,  Allan- 
factor,  detrended-fluctuation-analysis,  and  wavelet-transform  methods. 

Receiver-operating-characteristic  analysis  was  used  to  compare  and  contrast  the  per¬ 
formance  of  these  sixteen  measures  with  respect  to  their  abilities  to  properly  distinguish 
heart-failure  patients  from  normal  subjects  over  a  broad  range  of  record  lengths  (minutes  to 
hours).  Scale-dependent  measures  rendered  performance  that  was  substantially  superior  to 
that  of  scale-independent  ones.  Judged  on  the  basis  of  performance  and  computation  time, 
two  measures  of  the  sixteen  emerged  at  the  top  of  the  list:  the  wavelet-transform  standard 
deviation  crwav(32)  and  the  RR-interval-based  power  spectral  density  SV(l/32).  They  share 
in  common  the  dependence  on  a  single  scale,  or  small  range  of  scales,  near  32  heartbeat 
intervals.  The  behavior  of  the  ECG  at  this  particular  scale,  corresponding  to  about  25  sec, 
turns  out  to  be  a  significant  marker  for  the  presence  of  cardiac  dysfunction.  The  scale  value 
itself  is  far  more  important  than  the  details  of  the  measures  used  to  examine  it. 

Application  of  these  techniques  to  a  large  database  of  normal  and  CHF  records  will  make 
it  possible  to  uncover  just  how  the  ECGs  of  pathological  patients  differ  from  those  of  normal 
subjects  in  the  vicinity  of  this  special  scale.  This  would  facilitate  the  development  of  an  op¬ 
timal  diagnostic  measure,  which  might  take  the  form  of  a  specialized  weighting  function  over 
the  interval-based  power  spectral  density  or  over  some  combination  of  elementary  measures. 

We  also  addressed  an  issue  of  fundamental  importance  in  cardiac  physiology:  the  deter¬ 
mination  of  whether  the  RR  time  series  arises  from  a  chaotic  attractor  or  has  an  underlying 
stochastic  origin.  Using  nonlinear-dynamics  theory,  together  with  surrogate-data  analysis, 
we  established  that  the  RR  sequences  from  both  normal  subjects  and  pathological  patients 
have  stochastic,  rather  than  deterministic,  origins.  The  use  of  a  special  embedding  of  the 
differences  between  adjacent  RR  intervals  enabled  us  to  reach  this  conclusion.  This  tech¬ 
nique  substantially  reduced  the  natural  correlations  inherent  in  the  interbeat-interval  time 
series  which  can  masquerade  as  chaos  and  confound  simpler  analyses. 

Finally,  we  developed  a  jittered  integrate-and-fire  model  built  around  a  fractal-Gaussian- 
noise  kernel.  Simulations  based  on  this  model  provided  realistic,  though  not  perfect,  replicas 
of  real  heartbeat  sequences.  The  model  was  least  successful  in  predicting  interbeat-interval 
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statistics  and  fractal-exponent  values.  Nevertheless  it  could  find  use  in  some  applications 
such  as  pacemaker  excitation. 

To  confirm  the  observations  we  have  reported,  it  will  be  desirable  to  carry  out  continua¬ 
tion  studies  using  large  numbers  of  out-of-sample  records.  It  will  also  be  useful  to  examine 
how  the  performance  of  the  various  measures  correlates  with  the  severity  of  cardiac  dysfunc¬ 
tion  as  judged,  for  example,  by  the  NYHA  clinical  scale.  Comparisons  with  existing  CHF 
studies  that  make  use  of  HRV  measures  should  also  be  conducted. 
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Appendix  A 


Table  Al  provides  a  summary  of  some  significant  statistics  of  the  RR  records  analyzed  in 
this  Chapter,  before  truncation  to  75821  intervals. 
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Figure  Captions 


Fig.  1.  In  HRV  analysis  the  electrocardiogram  (ECG)  schematized  in  (a)  is  represented  by 
a  sequence  of  times  of  the  R  phases  which  form  an  unmarked  point  process  [vertical  arrows 
in  (b)].  This  sequence  may  be  analyzed  as  a  sequence  of  counts  (iVj}(T)  in  a  predetermined 
time  interval  T  as  shown  in  (c),  or  as  a  sequence  of  interbeat  (RR)  intervals  {t*}  as  shown  in 
(d).  The  sequence  of  counts  forms  a  discrete-time  random  counting  process  of  nonnegative 
integers  whereas  the  sequence  of  intervals  forms  a  sequence  of  positive  real-valued  random 
numbers. 

Fig.  2.  Estimating  the  wavelet  transform  using  the  Haar  wavelet:  (a)  original  Haar 
wavelet;  (b)  delayed  and  scaled  version  of  the  wavelet  (m  =  16,  n  =  3);  and  (c)  time  series 
multiplied  by  this  wavelet. 

Fig.  3.  (a)  Series  of  interbeat  intervals  r*  versus  interval  number  i  for  a  typical  nor¬ 
mal  patient  (data  set  nl6265).  (Adjacent  values  of  the  interbeat  interval  are  connected  by 
straight  lines  to  facilitate  viewing.)  Substantial  trends  are  evident.  The  interbeat-interval 
standard  deviation  crint  =  SDNN  is  indicated,  (b)  Wavelet  transform  Wm:n(m )  (calculated 
using  a  Daubechies  10-tap  analyzing  wavelet)  at  three  scales  (m  =  4,  16,  256)  vs.  interval 
number  i  ( —mn )  for  the  RR-interval  sequence  shown  in  (a).  The  trends  in  the  original 
interbeat-interval  time  series  are  removed  by  the  wavelet  transformation.  The  wavelet- 
transform  standard  deviation  crwav(m)  for  this  data  set  is  seen  to  increase  with  the  scale 
m. 


Fig.  4.  Haar-wavelet-transform  standard  deviation  <rwav(m)  vs  scale  m  for  the  12  normal 
subjects  (+),  12  CHF  patients  without  (s)  atrial  fibrillation  (x),  and  the  3  CHF  patients 
with  (c)  atrial  fibrillation  (A).  Each  data  set  comprises  the  first  75821  RR  intervals  of  a 
recording  drawn  from  the  Beth-Israel  Hospital  heart-failure  database.  Complete  separation 
of  the  normal  subjects  and  heart-failure  patients  is  achieved  at  scales  m  =  16  and  32  interbeat 
intervals. 

Fig.  5.  Individual  value  plots  (data)  for  the  16  measures.  Each  panel  corresponds  to  a 
different  HRV  measure.  The  seven  columns  in  each  panel,  from  left  to  right,  comprise  data 
for  (1)  12  normal  subjects  (+),  (2)  12  CHF  patients  s  AF  (x),  (3)  3  CHF  patients  c  AF  (A), 
(4)  1  heart-transplant  patient  (^),  (5)  1  ventricular  tachycardia  patient  (<0>),  (6)  1  sudden- 
cardiac-death  patient  (□),  and  (7)  2  sleep  apnea  patients  (+).  Each  data  set  comprises  75821 
RR  intervals  except  for  the  two  sleep  apnea  data  sets  which  comprise  16874  and  15751  RR 
intervals  respectively.  The  six  measures  highlighted  in  boldface  font  succeed  in  completely 
separating  normal  subjects  and  CHF  patients  (s  and  c  atrial  fibrillation)  VLF.  LF.  A(10), 
<rwav(32),  5  (1/32),  and  DFA(32). 

Fig.  6.  Positive  (solid  curves)  and  negative  (dotted  curves)  predictive  values  for  all 
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16  HRV  measures  plotted  against  the  threshold  6 ,  each  in  its  own  panel.  These  curves 
are  constructed  using  the  12  normal  and  12  heart-failure  (s  AF)  records  drawn  from  the 
CHF  database,  each  of  which  has  been  truncated  to  75821  RR  intervals.  The  six  measures 
highlighted  in  boldface  font  exhibit  threshold  regions  for  which  both  the  positive  and  negative 
predictive  values  are  unity:  VLF.  LF.  A(10),  <rwav(32),  S  (1/32),  and  DFA(32).  This 
indicates  that  the  normal  subjects  and  CHF  (s  AF)  patients  can  be  completely  distinguished 
by  these  six  measures,  in  accordance  with  the  results  established  in  Fig.  5. 

Fig.  7.  Receiver-operating-characteristic  (ROC)  curves  (sensitivity  vs  1  —  specificity)  for 
all  16  HRV  measures,  each  in  its  own  panel.  These  curves  are  constructed  using  the  12 
normal  and  12  heart-failure  (s  AF)  records  drawn  from  the  CHF  database,  each  of  which 
has  been  truncated  to  75821  RR  intervals.  The  six  measures  highlighted  in  boldface  font 
simultaneously  achieve  100%  sensitivity  and  100%  specificity  so  that  the  ROC  curve  is  per¬ 
fectly  square:  VLF.  LF.  A(10),  erwav(32),  S  (1/32),  and  DFA(32).  This  indicates  that 
the  normal  subjects  and  CHF  (s  AF)  patients  can  be  completely  distinguished  by  these  six 
measures,  in  accordance  with  the  results  established  in  Figs.  5  and  6. 

Fig.  8.  Diagnostic  accuracy  (area  under  ROC  curve)  as  a  function  of  record  length 
analyzed  L  (number  of  RR  intervals)  for  the  16  HRV  measures  (mean  ±1  S.D.).  An  area 
of  unity  corresponds  to  perfect  separability  of  the  two  classes  of  patients.  The  six  measures 
highlighted  in  boldface  font  (VLF,  LF.  A(10),  erwav(32),  S  (1/32),  and  DFA(32)  ) 
provide  such  perfect  separability  at  the  longest  record  lengths,  in  accordance  with  the  results 
in  Fig.  7.  As  the  record  length  decreases  performance  degrades  at  a  different  rate  for  each 
measure.  The  five  scale-independent  measures,  an,  am,  and  cur,  perform  poorly  at 

all  data  lengths. 

Fig.  9.  1— ROC  area  as  a  function  of  data  length  L,  on  a  doubly  logarithmic  plot, 

for  the  six  most  effective  measures:  VLF,  A(10),  DFA(32),  crwav(32),  AT(l/32),  and  LF. 
The  unconventional  form  of  this  ROC  plot  allows  small  deviations  from  unity  area  to  be 
highlighted.  The  trends  exhibited  by  all  six  measures  are  broadly  similar,  but  VLF,  A(10), 
and  DFA(32)  exhibit  more  degradation  at  shorter  data  lengths  than  do  the  other  three 
measures.  Thus  three  measures  emerge  as  superior:  LF.  crwav(32),  and  S  (1/32).  When 
computation  time  is  taken  into  account  in  accordance  with  the  results  provided  in  Table  2, 
<Twav(32)  and  S  (1/32)  are  the  two  preferred  measures. 

Fig.  10.  Capacity  dimension  Dq  as  a  function  of  embedding  dimension  p  for  the  12  normal 
subjects.  The  three  curves  shown  for  each  subject  correspond  to  the  original  differenced 
intervals  (solid  curves),  the  shuffled  surrogates  (dashed  curves),  and  the  phase-randomized 
surrogates  (dotted  curves).  The  dashed  curves  are  quite  similar  to  the  solid  curves  in  each 
panel  while  the  12  dotted  curves  closely  resemble  each  other. 

Fig.  11.  Capacity  dimension  Dq  as  a  function  of  embedding  dimension  p  for  the  12  CHF 
patients  without  atrial  fibrillation.  The  three  curves  shown  for  each  patient  correspond  to 
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the  original  differenced  intervals  (solid  curves),  the  shuffled  surrogates  (dashed  curves),  and 
the  phase-randomized  surrogates  (dotted  curves).  The  dashed  curves  are  reasonably  similar 
to  the  solid  curves  in  most  panels  while  the  12  dotted  curves  closely  resemble  each  other. 

Fig.  12.  Capacity  dimension  Dq  as  a  function  of  embedding  dimension  p  for  the  eight 
patients  with  other  cardiac  pathologies.  The  three  curves  shown  for  each  patient  correspond 
to  the  original  differenced  intervals  (solid  curves),  the  shuffled  surrogates  (dashed  curves), 
and  the  phase-randomized  surrogates  (dotted  curves).  The  dashed  curves  are  reasonably 
similar  to  the  solid  curves  in  most  panels  while  the  8  dotted  curves  closely  resemble  each 
other. 

Fig.  13.  Simple  integrate-and-fire  model  often  used  in  cardiology  and  neurophysiology. 
A  rate  function  A (t)  is  integrated  until  it  reaches  a  preset  threshold  6  whereupon  a  point 
event  is  generated  and  the  integrator  is  reset  to  zero. 

Fig.  14.  Schematic  representation  of  the  jittered  integrate-and-fire  model  for  generating 
a  simulated  RR-interval  series.  An  underlying  rate  process  A (t),  assumed  to  be  bandlimited 
fractal  Gaussian  noise,  is  integrated  until  it  reaches  a  fixed  threshold  9 ,  whereupon  a  point 
event  is  generated.  The  occurrence  time  of  the  point  event  is  jittered  in  accordance  with 
a  Gaussian  distribution  and  the  integrator  is  then  reset.  The  continuous  rate  process  is 
thereby  converted  into  a  point  process  representing  the  sequence  of  R  phases  in  the  human 
heartbeat. 

Fig.  15.  Comparison  between  data  (solid  curves)  and  simulations  (dotted  curves)  for  a 
single  normal  subject  (left  column,  data  set  nl6265),  and  a  single  CHF  s  AF  patient  (right 
column,  data  set  a6796).  The  parameters  E[rj]  (representing  the  mean  interbeat  interval) 
and  aw  (representing  the  fractal  exponent)  used  in  the  simuluations  were  derived  from  the 
data  for  the  two  individuals.  The  normal-subject  simulation  parameters  were  1/A  =  0.74  sec, 
as  =  1.0,  /o  =  0.0012  cycles/sec,  and  a  =  0.022;  the  CHF-patient  simulation  parameters 
were  1/A  =  0.99  sec,  as  =  1.5,  fo  =  0.00055  cycles/sec,  and  a  =  0.025.  The  lengths 
of  the  simulated  data  sets  were  identical  to  those  of  the  experimental  records  analyzed, 
(a)  RR-interval  sequence  over  the  entire  data  set.  Qualitative  agreement  is  apparent  in 
both  the  normal  and  heart-failure  panels,  (b)  Wavelet-transform  standard  deviation  versus 
scale.  The  simulations  mimic  the  scaling  properties  of  the  data  in  both  cases,  as  well  as  the 
gentle  flattening  of  erwav  for  the  heart-failure  patient  at  small  scale  values,  (c)  Interbeat- 
interval  histogram.  The  model  is  satisfactory  for  the  normal  subject  but  fails  to  capture  the 
narrowing  of  the  histogram  (reduction  of  crint)  for  the  heart-failure  patient,  (d)  Spectrum  of 
the  sequence  of  RR  intervals  (the  data  are  displaced  upward  by  a  factor  of  102  for  clarity). 
The  simulations  capture  the  subtleties  in  the  spectra  quite  well,  including  the  whitening  of 
the  heart-failure  spectrum  at  high  frequencies. 

Fig.  16.  Individual  value  plots  (jittered  integrate-and-fire  simulations)  for  the  16  mea¬ 
sures.  Each  panel  corresponds  to  a  different  HRV  measure.  The  three  columns  in  each  panel, 
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from  left  to  right,  comprise  data  for  (1)  12  simulations  using  normal-subject  parameters  (+), 
(2)  12  simulations  using  CHF-patient  (s  AF)  parameters  (x),  and  (3)  3  simulations  using 
CHF-patient  (c  AF)  parameters  (A).  Each  simulation  comprises  75821  RR  intervals  and  is 
carried  out  using  parameters  drawn  from  a  single  actual  data  set.  The  five  measures  high¬ 
lighted  in  boldface  font  succeed  in  completely  separating  normal-subject  and  CHF-patient  (s 
and  c  atrial  fibrillation)  simulations:  VLF.  LF.  <rwav(32),  S  (1/32),  and  DFA(32).  Each 
panel  should  be  compared  with  the  corresponding  panel  for  the  actual  normal  and  CHF  data 
in  Fig.  5  (leftmost  three  columns). 

Fig.  17.  Area  under  the  ROC  curve  (jittered  integrate-and-fire  simulations)  as  a  function 
of  number  of  RR  intervals  analyzed  (L)  for  the  16  HRV  measures  (mean  ±1  S.D.).  An  area  of 
unity  corresponds  to  perfect  separability  of  the  two  classes  of  simulations.  The  five  measures 
highlighted  in  boldface  font  [VLF.  LF.  <xwav(32),  S  (1/32),  and  DFA(32)]  provide  such 
perfect  separability,  but  only  for  the  largest  number  of  RR  intervals  analyzed.  Each  panel 
should  be  compared  with  the  corresponding  panel  for  the  actual  normal  and  CHF  data  in 
Fig.  8.  The  ROC-area  simulations  differ  from  those  for  the  data  in  two  principal  respects: 
the  performance  of  the  two  interval  measures  SDANN  and  <7int  severely  degrades  as  the 
number  of  RR  intervals  increases,  and  the  performance  of  the  five  fractal-exponent  measures 
is  substantially  enhanced  as  the  number  of  RR  intervals  increases. 

Fig.  18.  Simulation  accuracy  for  four  measures  with  their  correlation  coefficients  p.  The 
jittered  integrate-and-fire  model  performs  remarkably  well  for  the  measures  <SV(l/32)  and 
c wav (32);  however  it  does  not  perform  nearly  as  well  for  the  measures  crint  and  aw-  These 
disagreements  highlight  problems  with  the  simple  jittered  integrate-and-fire  model. 

Fig.  19.  Potential  modifications  to  the  simple  integrate-and-fire  model.  Multiple  phys¬ 
iologically  based  inputs  and  adaptive  feedback  could  serve  to  produce  more  realistic  RR- 
interval  simulations. 
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Tables 


TABLE  1:  Detection  distances  h  and  d  for  all  16  HRV  measures  applied  to  the  12  normal 
and  12  CHF  records  of  length  Lmax.  The  measures  are  listed  in  order  of  descending  value 
of  h.  The  rankings  according  to  d  differ  from  those  according  to  h.  The  wavelet-transform 
standard  deviation  <rwav(32)  and  variance  <r^av(32),  though  related  by  a  simple  monotonic 
transformation,  yield  different  values  of  h  and  have  different  rankings. 


Measure  h  d 


DFA(32) 
^"wav  (32) 

-4(10) 

VLF 

St(1/32) 

C^int 

ffj„(32) 

OiD 

SDANN 

LF 

pNN50 

LF/HF 

aw 

OCR 

HF 

as 

OtA 


2.48253 

2.33614 

2.32522 

1.84285 

1.77422 

1.74750 

1.71343 

1.64883 

1.46943 

1.36580 

1.36476 

1.24507 

1.09916 

1.02367 

0.85361 

0.82125 

0.38778 


1.81831 

1.70153 

1.77482 

1.56551 

1.55200 

1.32475 

1.47165 

1.17679 

1.04079 

1.28686 

1.20896 

0.91444 

0.77800 

0.72463 

0.73077 

0.58071 

0.27895 
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TABLE  2:  Computation  times  (to  the  nearest  10  msec)  for  the  16  HRV  measures  for  data 
sets  comprising  75821  RR  intervals.  The  long  execution  times  for  the  two  DFA  measures 
results  from  the  fact  that  it  is  an  N2  process  whereas  the  14  other  methods  are  either  N  or 
iVlog(iV). 


Execution 

Measure  Time  (msec) 


VLF,  LF,  HF,  and  LF/HF 

330 

pNN50 

40 

SDANN 

160 

C^int 

190 

-4(10) 

160 

dwav(32) 

20 

<SV(l/32) 

60 

DFA(32) 

650,090 

OLD 

650,110 

Olw 

220 

as 

920 

OLA 

610 

olr 

570 
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TABLE  Al:  Elementary  statistics  of  the  original  RR-interval  recordings,  before  truncation 
to  75821  RR  intervals. 


File¬ 

name 

Condition 

Number  of 
RR  Intervals 

Recording 
Dnr.  (sec) 

RR-Interval 
Mean  (sec) 

Mean  Rate 
(sec)-1 

RR-interval 
S.D.  (sec) 

nl6265 

Normal 

100460 

80061 

0.797 

1.255 

0.171 

nl6272 

93177 

84396 

0.906 

1.104 

0.142 

nl6273 

')') 

89846 

74348 

0.828 

1.208 

0.146 

nl6420 

V 

102081 

77761 

0.762 

1.313 

0.101 

nl6483 

')') 

104338 

76099 

0.729 

1.371 

0.089 

nl6539 

5) 

108331 

84669 

0.782 

1.279 

0.150 

nl6773 

99 

82160 

78141 

0.951 

1.051 

0.245 

nl6786 

99 

101630 

84052 

0.827 

1.209 

0.116 

nl6795 

99 

87061 

74735 

0.858 

1.165 

0.212 

nl7052 

99 

87548 

76400 

0.873 

1.146 

0.159 

nl7453 

99 

100674 

74482 

0.740 

1.352 

0.103 

nc4 

99 

88140 

71396 

0.810 

1.235 

0.132 

a6796 

HF  s  AF 

75821 

71934 

0.949 

1.054 

0.091 

a8519 

99 

80878 

71948 

0.890 

1.124 

0.062 

a8679 

99 

119094 

71180 

0.598 

1.673 

0.051 

a9049 

99 

92497 

71959 

0.778 

1.285 

0.058 

a9377 

99 

90644 

71952 

0.794 

1.260 

0.060 

a9435 

99 

114959 

71158 

0.619 

1.616 

0.034 

a9643 

99 

148111 

71958 

0.486 

2.058 

0.024 

a9674 

99 

115542 

71968 

0.623 

1.605 

0.084 

a9706 

99 

115064 

71339 

0.620 

1.613 

0.100 

a9723 

99 

115597 

71956 

0.622 

1.607 

0.027 

a9778 

99 

93607 

71955 

0.769 

1.301 

0.070 

a9837 

99 

115205 

71944 

0.624 

1.601 

0.066 

a7257 

HF  c  AF 

118376 

71140 

0.601 

1.664 

0.039 

a8552 

99 

111826 

71833 

0.642 

1.557 

0.062 

a8988 

99 

118058 

71131 

0.603 

1.660 

0.091 

tp987  TRANSPLANT  106394 

67217 

0.632 

1.583 

0.083 

vt973 

VT 

86992 

70044 

0.805 

1.242 

0.202 

sd984 

SCD 

77511 

62786 

0.810 

1.234 

0.089 

slp59 

APNEA 

16874 

14399 

0.853 

1.172 

0.103 

slp66 

99 

15751 

13199 

0.838 

1.193 

0.103 
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